{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "502_Convergent.ipynb", "provenance": [], "toc_visible": true, "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "gbK7tnuaP7Cc" }, "source": [ "# Convergence of a Multistep Method\n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "## Overview\n", "A one-step or multistep method is used to approximate the solution of an initial value problem of the form\n", "\\begin{equation} \\frac{dy}{dt}=f(t,y),\\end{equation}\n", "with the initial condition\n", "\\begin{equation} y(a)=\\alpha.\\end{equation}\n", "The method should only be used if it satisfies the three criteria:\n", "1. that difference equation is __consistent__ with the differential equation;\n", "2. that the numerical solution __convergent__ to the exact answer of the differential equation;\n", "3. that the numerical solution is __stable__.\n", "\n", "In the notebooks in this folder we will illustate examples of consistent and inconsistent, convergent and non-convergent, and stable and unstable methods. \n" ] }, { "cell_type": "markdown", "metadata": { "id": "dT1-22QhP7Cd" }, "source": [ "## Introduction to Convergence\n", "\n", "In this notebook we will illustate an non-convergent method. The video below outlines the notebook.\n" ] }, { "cell_type": "code", "metadata": { "id": "RFMbOXWdP7Cf", "outputId": "0085863c-8389-46cc-de8e-67532a89b00c", "colab": { "base_uri": "https://localhost:8080/", "height": 336 } }, "source": [ "from IPython.display import HTML\n", "HTML('')" ], "execution_count": 1, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "yBuEvQltP7Ck" }, "source": [ "### Definition\n", "\n", "The solution of a multi-step methods $w_i$ is said to be __convergent__ with respect to\n", "the exact solution of the differential equation if\n", "\\begin{equation} \\max_{h \\rightarrow 0}\\max_{1 \\leq i \\leq N}|y(t_i)-w_i|=0.\\end{equation}\n", "All the Runge Kutta, Adams-Bashforth and Adams-Moulton methods are convergent." ] }, { "cell_type": "markdown", "metadata": { "id": "r8ft_ye7P7Cl" }, "source": [ "## 2-step Abysmal Butler Multistep Method \n", "\n", "This notebook will illustrate a non-convergent multistep method using the Abysmal-Butler method, named with great pride after the author.\n", "The 2-step Abysmal Butler difference equation is given by\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(4f(t_i,w_i)-3f(t_{i-1},w_{i-1})),\\end{equation}\n", "The final section of this notebooks shows that the method is non-convergent for all differential equations.\n", "\n", "## Intial Value Problem\n", "To illustrate convergence we will apply Abysmal-Butler multistep method to the linear intial value problem\n", "\\begin{equation} y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2) \\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1,\\end{equation}\n", "with the exact solution\n", "\\begin{equation}y(t)= 2e^{-t}+t-1.\\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "id": "xCeYXckcP7Cm" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "jZ0cRNZVP7Cn" }, "source": [ "import numpy as np\n", "import math \n", "import pandas as pd\n", "\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")" ], "execution_count": 2, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Dq46Ycx6P7Cq" }, "source": [ "### Defining the function\n", "\\begin{equation}f(t,y)=t-y.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "vs6cXGsKP7Cr" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": 3, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "j9Gf-5FVP7Cu" }, "source": [ "## Discrete Interval\n", "To illustrtate that the method is internally convergent but not convergent with the exact solution we define two discrete meshes, one coarse and one fine.\n", "### Coarse mesh\n", "Defining the step size $h$ from the interval range $a \\leq t \\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b - a}{N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "$$t_i=t_0+ih,$$\n", "where $t_0=a,$ for $i=0,1...,N$.\n", "### Fine mesh\n", "Defining the step size $h/2$ from the interval range $a≤t≤b$ and number of steps $2N$ \n", "\\begin{equation}h=\\frac{b−a}{2N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation}\n", "where $t_0=a,$ for $i =0,1,...2N$.\n", "\n", "The plot below shows the coarse (red) and fine (green) discrete time intervals over the domain." ] }, { "cell_type": "code", "metadata": { "id": "64MsysByP7Cu", "outputId": "70d46f79-4b1b-4ba8-8f29-655b51c6c973", "colab": { "base_uri": "https://localhost:8080/", "height": 281 } }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=16\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "N2=N*2\n", "h2=(b-a)/(N2)\n", "t2=np.arange(a,b+h2,h2)\n", "w2=np.zeros(len(t2))\n", "\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t2,0.01*np.ones(len(t2)),'o:',color='green',label='Fine Mesh')\n", "plt.plot(t,0*t,'o:',color='red',label='Coarse Mesh')\n", "plt.xlim((0,2))\n", "plt.ylim((-0.1,.1))\n", "\n", "plt.legend()\n", "plt.title('Illustration of discrete time points')\n", "plt.show()" ], "execution_count": 4, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAEICAYAAAA0vXKRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xVZfn38c9XQAgPiIpKoqDlIfCUTKaVOuUJTcVSCx80NJU8pPFL7cGoNJXUTpJP9TMyxZQEo0yoFAEl88yMIQePgAcg0OHgKCIIej1/rDXjZth7DuwNM2vzfb9e+zX7vte97n2te9+z55p12EsRgZmZmZll0xatHYCZmZmZbTgnc2ZmZmYZ5mTOzMzMLMOczJmZmZllmJM5MzMzswxzMmdmZmaWYU7mzNowSWdLejSnHJI+2ZoxFSLpFkk/bIXXvVDSG5JWSNqhGe1flXR0+vz7km7d+FFuPJIGSnqwteNojKTd0/enXWvHYlaOnMyZtTG5ycZG6n+UpOuK7GOdJBMgIi6IiGuLi67FcXQAfgkcGxFbR8TSlqwfET+JiPM2TnTrk9QrTcjbl2r9iBgdEceWLsrSi4jX0/fng6baFjtGZpsjJ3Nmto6M/RHdGegEzG7tQCBzY2dmZcLJnFlGSZoq6byccv3eMiVukvSmpLclzZS0n6TBwEDge+lhrwlp+1cl/V9JM4B3JbWXNFTSXEnvSHpO0lfStp8CbgEOS/t4K61fZ4+fpPMlzZG0TNJ4SR/PWRaSLpD0sqS3JP1GkgpsZ0dJIyT9N32MSOv2Bl5Mm70l6aEC658l6TVJSyUNa7Dsakl3pc87SborbfeWpGmSdk6XbS/p9vT1l0v6W1pfKWlBOnaLgdslbZEzdksl3SNp+/QlH8mJd4Wkw9J+vinp+bTviZJ6Fnjb11u/wKH4i9KxfUfStZI+IenxdC7cI2nLnPYnSpqebvPjkg4o8Np1fV8qaZ6kJZJ+JmmLdNkWkn6QjvWbkv4oqUu6bJ29bencvVbSY2mMD0rasZFt/KSkf0mqTV93bKEYzTZHTubMytOxwBHA3kAX4GvA0ogYCYwGfpoe9jopZ50zgC8D20XEWmAucHi6/o+BuyR1j4jngQuAJ9I+tmv44pK+BFyfvm534DVgTINmJwKfAQ5I2x1XYFuGAYcCBwEHAocAP4iIl4A+aZvtIuJLeeLoDfwvcBbwcWAHoEeB1xmUbutuabsLgPfSZXcCndPX2wm4KWe9XYDtgZ7AYOAS4BTgyPQ1lwO/SdsekRPv1hHxhKT+wPeBrwLdgH8DdxeIcb31C7Q7DuhLMm7fA0YCZ6bbth/Je42kTwO3Ad9Kt/l3wHhJHQv0C/AVoAI4GOgPfDOtPzt9fBHYE9ga+HUj/fwf4ByS8dwSuLyRbbwWeBDoSvL+/b9G+jXb7DiZMytPa4BtgH0BRcTzEbGoiXVujoj5EfEeQET8OSL+GxEfRsRY4GWSRKo5BgK3RcQzEbEauJJkT16vnDY3RMRbEfE68DBJslaor2si4s2IqCFJLM9qZhynAX+PiEfSOH4IfFig7RqShOaTEfFBRFRHxNuSugPHAxdExPKIWBMR/8pZ70PgqohYnY7dBcCwiFiQvubVwGkqfAj2AuD69D1aC/wEOKiRvXPN8dOIeDsiZgOzgAcjYl5E1AL3A59O2w0GfhcRT6XbfAewmiQJLOTGiFiWvm8jSBNDkvfpl+nrrCB5zwc0st23R8RL6ZjdQ+H3H5L3pifw8YhYFRGPNtLWbLPjZM6sDEXEQyR7RX4DvClppKRtm1htfm5B0jdyDr+9RbJHZ8f8q67n4yR74+riWQEsBXbNabM45/lKkj05TfaVPv94gbb51q3froh4N40jnzuBicCY9HDqT5VcYLEbsCwilhdYryYiVuWUewL35ozb88AHJOf35dMT+FVO+2WAWHesWuqNnOfv5SnXjXVP4LK6105ffzcaH9/ceZL7XuR7n9pTeLub+/5DsndRwNOSZkv6ZiNtzTY7TubMsutdkkN/dXbJXRgRN0dEX6A3yeHWK+oWFeivvj7dK/R74NvADumh1Fkkf1Ab66POf0kShbr+tiLZ67WwifWa7AvYPa1rjkUkyUldHJ3TONaT7nH7cUT0Bj5Hchj4GyTJy/aS1jucXLdqg/J84PiI2C7n0SkiFuZpW9f+Ww3afywiHm/GaxVrPjC8wWt3johCh3khZzxZ973I9z6tZd1EsjnW28aIWBwR50fEx0kOCf9WbfQresxag5M5s+yaDnxVUuf0D9u5dQskfUbSZ9M9S+8Cq/jo8OIbJOc0NWYrkj+qNWl/55DsmavzBtAj90T6Bu4GzpF0UHr+1U+ApyLi1ZZsYE5fP5DULT1J/kfAXc1cdxxwoqQvpLFeQ4HPPUlflLS/ku9Ce5vk0N6H6eHp+0kSiK6SOkg6Il8fqVuA4XWHSdO4+6fLakjehz0btL9SUp+0fRdJpxfoO9/6xfg9cEE6VyRpK0lflrRNI+tckY7DbsB3gLqLEe4G/kfSHpK2JnnPx6aHjltivW2UdLqkunMdl5PMzUKHy802O07mzLLrJuB9ksTqDpILG+psS/KHejnJ4a6lwM/SZX8AeqeH1f6Wr+OIeA74BfBE2v/+wGM5TR4i+TqQxZKW5Fl/Msn5aX8h2Tv2CWDABm0lXAdUATOAmcAzaV2T0nPGLgb+lMaxHFhQoPkuJMnf2ySHRv9FcugVknP01gAvAG8CQxp52V8B44EHJb0DPAl8No1nJTAceCwd/0Mj4l7gRpLDu2+T7AE9vsD2rLd+k4PQiIioAs4nOSS/HJhDchFDY+4Dqkn+mfgHyXyC5EKKO0muRn2F5B+ISzYgpnzb+BngKUkrSMb2OxExr6V9m5UrRZR6r72ZmZUjSQHsFRFzWjsWM/uI98yZmZmZZVhJkjlJ/SS9qOQLQofmWX6EpGckrZV0WoNlg5R8ueXLkgbl1PdV8kWncyTdLOX/QlEzMzOzzVnRyVx6svBvSM7x6A2ckX5RZ67XSc7D+FODdbcHriI5n+QQ4CpJXdPF/0tyLsde6aNfsbGamdmGiwj5EKtZ21OKPXOHAHPSL4p8n+Rb3vvnNoiIVyNiButffXQcMCn9AsrlwCSgX/olndtGxJORnNT3R5JvVDczMzOzHKW4KfSurPslkgtIr9zawHV3TR8L8tSvR8m9JgcDbLXVVn333XffZr60mZmZWeuprq5eEhHdiu2nFMlcq0rvNTkSoKKiIqqqqlo5IjMzM7OmSXqt6VZNK8Vh1oWs+43gPWj+t7wXWnch694MuyV9mpmZmW02SpHMTQP2Sr/1e0uSLwYd38x1JwLHpt8m3hU4FpiYfuP625IOTa9i/QbJF1WamZmZWY6ik7n0Vi3fJknMngfuiYjZkq6RdDLU31poAXA68DtJs9N1lwHXkiSE04Br0jqAi4BbSb6RfC7J7XTMzMzMLEdZ3QHC58yZmZkl1qxZw4IFC1i1alVrh7LZ69SpEz169KBDhw7r1EuqjoiKYvvP/AUQZmZmtr4FCxawzTbb0KtXL/y9+60nIli6dCkLFixgjz322Civ4dt5mZmZlaFVq1axww47OJFrZZLYYYcdNuoeUidzZmZmZcqJXNuwsd8HJ3NmZmZmGeZkzszMzDaKdu3acdBBB9U/Xn31VT73uc+VpO+rr74aScyZ89HtgkeMGIEkNuRiyLPPPptx48aVJLZNzcmcmZmZMXrmaHqN6MUWP96CXiN6MXrm6KL7/NjHPsb06dPrH7169eLxxx8vQbSJ/fffnzFjxtSX//znP9OnT5+S9Z8VTubMzMw2c6NnjmbwhMG8VvsaQfBa7WsMnjC4JAldQ1tvvTUAU6dOpbKyktNOO419992XgQMHUvd1adXV1Rx55JH07duX4447jkWLFuXt65RTTuG++5J7CsydO5cuXbqw44471i9/8MEHOeywwzj44IM5/fTTWbFiBQBDhw6ld+/eHHDAAVx++eX17R955BE+97nPseeee2ZqL52TOTMzs81A5ahKRk0fBcCaD9ZQOaqSu2bcBcCVk69k5ZqV67RfuWYlQx4YAsCSlUuoHFXJhBcnALB4xeJmveZ7771Xf4j1K1/5ynrL//Of/zBixAiee+455s2bx2OPPcaaNWu45JJLGDduHNXV1Xzzm99k2LBhefvfdttt2W233Zg1axZjxozh61//ev2yJUuWcN111zF58mSeeeYZKioq+OUvf8nSpUu59957mT17NjNmzOAHP/hB/TqLFi3i0Ucf5e9//ztDhw5t1ja2Bf6eOTMzs83cgrcX5K1funJpUf3WHWYt5JBDDqFHj+RW7HXn1G233XbMmjWLY445BoAPPviA7t27F+xjwIABjBkzhokTJzJlyhRuv/12AJ588kmee+45Pv/5zwPw/vvvc9hhh9GlSxc6derEueeey4knnsiJJ55Y39cpp5zCFltsQe/evXnjjTeK2vZNycmcmZnZZmDq2VPrn3do12Gd8u5ddue12tfWW2f3LrsDsGPnHddpv8vWu5Qkpo4dO9Y/b9euHWvXriUi6NOnD0888USz+jjxxBO54oorqKioYNttt62vjwiOOeYY7r777vXWefrpp5kyZQrjxo3j17/+NQ899NB68WTpDlk+zGpmZraZG37UcDp36LxOXecOnRl+1PBNHss+++xDTU1NfTK3Zs0aZs+eXbB9586dufHGG9c7FHvooYfy2GOP1V/t+u677/LSSy+xYsUKamtrOeGEE7jpppt49tlnN97GbCLeM2dmZraZG7j/QACGTRnG67Wvs3uX3Rl+1PD6+k1pyy23ZNy4cVx66aXU1taydu1ahgwZ0uhVqgMGDFivrlu3bowaNYozzjiD1atXA3DdddexzTbb0L9/f1atWkVE8Mtf/nKjbcumoiztRmxKRUVFbMh3y5iZmZWb559/nk996lOtHYal8r0fkqojoqLYvn2Y1czMzCzDnMyZmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyrCTJnKR+kl6UNEfSeve/kNRR0th0+VOSeqX1AyVNz3l8KOmgdNnUtM+6ZTuVIlYzMzOzclJ0MiepHfAb4HigN3CGpN4Nmp0LLI+ITwI3ATcCRMToiDgoIg4CzgJeiYjc+34MrFseEW8WG6uZmZltOosXL2bAgAF84hOfoG/fvpxwwgm89NJLrRrTqFGjkMTkyZPr6/72t78hiXHjxrW4v6uvvpqf//znpQyxxUqxZ+4QYE5EzIuI94ExQP8GbfoDd6TPxwFHSVKDNmek65qZmdmmNno09OoFW2yR/Bw9uqjuIoKvfOUrVFZWMnfuXKqrq7n++uuLvudpRPDhhx8W1cf+++/PmDEfpRx33303Bx54YFF9tqZSJHO7AvNzygvSurxtImItUAvs0KDN14GGN1C7PT3E+sM8yZ+ZmZmVwujRMHgwvPYaRCQ/Bw8uKqF7+OGH6dChAxdccEF93YEHHsjhhx9ORHDFFVew3377sf/++zN27FgAVqxYwVFHHcXBBx/M/vvvz3333QfAq6++yj777MM3vvEN9ttvP+bPn8/ZZ59dv/5NN90EwNy5c+nXrx99+/bl8MMP54UXXsgb2+GHH87TTz/NmjVrWLFiBXPmzOGggw6qX15dXc2RRx5J3759Oe6441i0aBEAN998M7179+aAAw5Y564Tzz33HJWVley5557cfPPNGzxmG6pN3M5L0meBlRExK6d6YEQslLQN8BeSw7B/zLPuYGAwwO67774pwjUzM8ueyko4++zksWYNHHMMnHcenHkmXHklrFy5bvuVK2HIEBg4EJYsgdNOg8sug5NOgsWLYZddGn25WbNm0bdv37zL/vrXvzJ9+nSeffZZlixZwmc+8xmOOOIIunXrxr333su2227LkiVLOPTQQzn55JMBePnll7njjjs49NBDqa6uZuHChcyalaQNb731FgCDBw/mlltuYa+99uKpp57ioosu4qGHHlrv9SVx9NFHM3HiRGprazn55JN55ZVXgOResJdccgn33Xcf3bp1Y+zYsQwbNozbbruNG264gVdeeYWOHTvWvybACy+8wMMPP8w777zDPvvsw4UXXkiHDh2a8aaURimSuYXAbjnlHmldvjYLJLUHugBLc5YPoMFeuYhYmP58R9KfSA7nrpfMRcRIYCQkt/MqakvMzMw2RwsW5K9fujR/fZEeffRRzjjjDNq1a8fOO+/MkUceybRp0zj++OP5/ve/zyOPPMIWW2zBwoUL6w/L9uzZk0MPPRSAPffck3nz5nHJJZfw5S9/mWOPPZYVK1bw+OOPc/rpp9e/Tt09WfMZMGAAN998M7W1tfziF7/gJz/5CQAvvvgis2bN4phjjgHggw8+oHv37gAccMABDBw4kFNOOYVTTjmlvq8vf/nLdOzYkY4dO7LTTjvxxhtv0KNHj9IOWiNKkcxNA/aStAdJ0jYA+D8N2owHBgFPAKcBD0V6U1hJWwBfAw6va5wmfNtFxBJJHYATgcmYmZnZhpk69aPnHTqsW9599+TQakN1R7x23HHd9k3slQPo06dPiy8oGD16NDU1NVRXV9OhQwd69erFqlWrANhqq63q23Xt2pVnn32WiRMncsstt3DPPfcwYsQItttuO6ZPn16o+3UccsghzJw5k86dO7P33nvX10cEffr04YknnlhvnX/84x888sgjTJgwgeHDhzNz5kwAOnbsWN+mXbt2rF27tkXbXayiz5lLz4H7NjAReB64JyJmS7pG0slpsz8AO0iaA3wXyP36kiOA+RExL6euIzBR0gxgOkmS+PtiYzUzM7M8hg+Hzp3XrevcOanfQF/60pdYvXo1I0eOrK+bMWMG//73vzn88MMZO3YsH3zwATU1NTzyyCMccsgh1NbWstNOO9GhQwcefvhhXsuXYAJLlizhww8/5NRTT+W6667jmWeeYdttt2WPPfbgz3/+M5AkZc8++2yjMd5www31e+Tq7LPPPtTU1NQnc2vWrGH27Nl8+OGHzJ8/ny9+8YvceOON1NbWsmLFig0en1IqyTlzEfFP4J8N6n6U83wVcHrD9dJlU4FDG9S9C+Q/0G5mZmalNXBg8nPYMHj99WSP3PDhH9VvAEnce++9DBkyhBtvvJFOnTrRq1cvRowYwRe+8AWeeOIJDjzwQCTx05/+lF122YWBAwdy0kknsf/++1NRUcG+++6bt++FCxdyzjnn1F/Vev311wPJnr0LL7yQ6667jjVr1jBgwIBGr1I9/vjj16vbcsstGTduHJdeeim1tbWsXbuWIUOGsPfee3PmmWdSW1tLRHDppZey3XbbbfD4lJLSo51loaKiIqqqqlo7DDMzs1b3/PPP86lPfaq1w7BUvvdDUnVEVBTbt2/nZWZmZpZhTubMzMzMMszJnJmZWZkqp1Opsmxjvw9O5szMzMpQp06dWLp0qRO6VhYRLF26lE6dOm2012gTd4AwMzOz0urRowcLFiygpqamtUPZ7HXq1GmjfomwkzkzM7My1KFDB/bYY4/WDsM2AR9mNTMzM8swJ3NmZmZmGeZkzszMzCzDnMyZmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLMCdzZmZmZhnmZM7MzMwsw5zMmZmZmWVYSZI5Sf0kvShpjqSheZZ3lDQ2Xf6UpF5pfS9J70manj5uyVmnr6SZ6To3S1IpYjUzMzMrJ0Unc5LaAb8Bjgd6A2dI6t2g2bnA8oj4JHATcGPOsrkRcVD6uCCn/n+B84G90ke/YmM1MzMzKzel2DN3CDAnIuZFxPvAGKB/gzb9gTvS5+OAoxrb0yapO7BtRDwZEQH8ETilBLGamZmZlZVSJHO7AvNzygvSurxtImItUAvskC7bQ9J/JP1L0uE57Rc00ScAkgZLqpJUVVNTU9yWmJmZmWVMa18AsQjYPSI+DXwX+JOkbVvSQUSMjIiKiKjo1q3bRgnSzMzMrK0qRTK3ENgtp9wjrcvbRlJ7oAuwNCJWR8RSgIioBuYCe6ftezTRp5mZmdlmrxTJ3DRgL0l7SNoSGACMb9BmPDAofX4a8FBEhKRu6QUUSNqT5EKHeRGxCHhb0qHpuXXfAO4rQaxmZmZmZaV9sR1ExFpJ3wYmAu2A2yJitqRrgKqIGA/8AbhT0hxgGUnCB3AEcI2kNcCHwAURsSxddhEwCvgYcH/6MDMzM7McSi4WLQ8VFRVRVVXV2mGYmZmZNUlSdURUFNtPa18AYWZmZmZFcDJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLMCdzZmZmZhnmZM7MzMwsw5zMmZmZmWWYkzkzMzOzDHMyZ2ZmZpZhTubMzMzMMszJnJmZmVmGOZkzMzMzyzAnc2ZmZmYZ5mTOzMzMLMOczJmZmZllmJM5MzMzswxzMmdmZmaWYSVJ5iT1k/SipDmShuZZ3lHS2HT5U5J6pfXHSKqWNDP9+aWcdaamfU5PHzuVIlYzMzOzctK+2A4ktQN+AxwDLACmSRofEc/lNDsXWB4Rn5Q0ALgR+DqwBDgpIv4raT9gIrBrznoDI6Kq2BjNzMzMylUp9swdAsyJiHkR8T4wBujfoE1/4I70+TjgKEmKiP9ExH/T+tnAxyR1LEFMZmZmZpuFUiRzuwLzc8oLWHfv2jptImItUAvs0KDNqcAzEbE6p+729BDrDyUp34tLGiypSlJVTU1NMdthZmZmljlt4gIISX1IDr1+K6d6YETsDxyePs7Kt25EjIyIioio6Nat28YP1szMzKwNKUUytxDYLafcI63L20ZSe6ALsDQt9wDuBb4REXPrVoiIhenPd4A/kRzONTMzM7McpUjmpgF7SdpD0pbAAGB8gzbjgUHp89OAhyIiJG0H/AMYGhGP1TWW1F7SjunzDsCJwKwSxGpmZmZWVopO5tJz4L5NciXq88A9ETFb0jWSTk6b/QHYQdIc4LtA3deXfBv4JPCjBl9B0hGYKGkGMJ1kz97vi43VzMzMrNwoIlo7hpKpqKiIqip/k4mZmZm1fZKqI6Ki2H7axAUQZmZmZrZhnMyZmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLsLJK5qr/W02vEb0YPXP0Bq0/euZoeo3oxRY/3mKD+2krfTiWth9LuW2PY2n7sZTb9jiWjdeHY9k0sdCdvhsURAPtrr766lL00yb8+Bc/vrp2v1oemPMAvbbrxQE7H9DsdUfPHM3gCYNZsnIJALWrW95PW+nDsbT9WMptexxL24+l3LbHsWwe21P2sVTD1Zdd/eNmB1BAWe2Zq7NyzUq+P+X7VI6q5K4Zd9XXVY6qZOyssQDUrqqlclQlf33+rwAMnTyUlWtWrtfPdx/4LgDza+dTOaqSyfMmAzBv+TwqR1Xyr1f/BcCLS17k3PvOzdvHsCnDmL54OpWjKpm+eDoA0xZOo3JUJbPeTG45+/j8x6kcVcn3Hvxe/jgmfpfKUZXMr50PwANzHqByVCWLVywGYMKLE6gcVVk/wYY8MCRvP1dOvhKAu2bcReWoStZ8sAaAUdNHUTmqsr7t76t/z3n3nVdwe3715K84+e6T6+t//vjPOfWeU+vLNzx6AwPGDQBg2JRhefv51oRv1ZevnHwlgycMri9f/uDlXPyPi+vLF/79woKxAAyeMLh+2wDOue8cfvTwj+rLZ/71TK7917UFYzl//Pn8/PGf19edfPfJ/OrJX9WXjx99PL+d9ttGt+c7938HgDUfrGly7jU2VxavWEzlqEoemPMA0PjcKxTLufedy7SF0wCanHuXTbwsbx/fm/Q9ACbPm9ysuVfod2jI/UOoHFVZv6yxuVdoe84bf159uTlzb/CEwQXH90cP/4hz7junvr7Q3CsUywV/v4DLH7y8vq6puddYLACn3nNqk3Pvkn9eUrCPylGVjJo+Cmh87jU2Vya8OAGgWXOvsbk7681ZVI6qbHLuXfHgFXn7uGziZVSOqmTe8nlA03Pvfx74n0Y/58bOGtvk3Gtse3477bccP/r4+vrG5l6h8c2dW03NvaY+5y7+x8VNzr2L/3FxwT4GjBvADY/eUF/f2NwrtD2X3n9pfbmpudfY2C5ZuYTKUZXNmnuNzd3H5z8O0OTcu/zBy/N/zj2YfM7969V/NWvuDZ1U4HPugeRzrnZVLdD43Mu3PcUqy2QOqH8Dmmvh2wvz1tesrGl2H6s/WJ23/vXa15vdx6IVi/LH8W7z4wBYunJp3voFby9odh+rPliVt74l29NY+3fXvNvsPt55/52NGst7a98ruo+l7+Uf83xKMVcaa1+o/3zefPfNvPWL3sk/Fwsp9DvUknEptD2r1uafi4UU+qBsyfgWarvi/RWbPJblq5YX3Ucp5kpj7VsSS10y1lChuVhI3T+vDbXkc25j/y625I92KT7nalfXFt1HY+2Xvbes2X20pc+5N1a8kbe+0N/cQha+U+BzrsDf3Hxauv3NUVb3ZtXHFaQ7e3p26cmrQ15t9rq9RvTitdrX1qtvST9tpQ/H0vZjKbftcSxtP5Zy2x7HsvH6cCybMJbfQfw31OwACijLPXOdO3Rm+FHDW7TO8KOG07lD56L6aSt9OJa2H0u5bY9jafuxlNv2OJaN14dj2fSxFC0iyuZBd6LnTT3jrhl3xYa4a8Zd0fOmnqGrtcH9tJU+HEvbj6XctsextP1Yym17HMvG68OxbJpY6E5ECfKfsjrMWlFREVVVVa0dhpmZmVmTJFVHREWx/ZTlYVYzMzOzzYWTOTMzM7MMczJnZmZmlmElSeYk9ZP0oqQ5kobmWd5R0th0+VOSeuUsuzKtf1HScc3tM6/qaujVC0Zv2C06SmL06CSGLbZwLI4lG3E4FsfiWMorlrYSh2NpMpa+lOZ2XsVfQQrtgLnAnsCWwLNA7wZtLgJuSZ8PAMamz3un7TsCe6T9tGtOn/kefSECIjp3jrhrw65SKcpddyWvXReHY3EsbT0Ox+JYHEt5xdJW4nAszYqlL6W5mrUUydxhwMSc8pXAlQ3aTAQOS5+3B5YAati2rl1z+sz36Jv7Ju2+e8SRR0bceWcyeLAxKmAAABJ2SURBVO++m5THjEnKb72VlP/yl6RcU5OUx49PyosWJeX770/Kr7+elCdNSspz5yblqVOT8gsvRHTsuO5EqXv07Bnxn/8k7f/zn6T9008n5Zkzk/JjjyXlF15IylOnJuW5c5PypElJ+fXXk/L99yflRYuS8vjxSbmmJinvuGP+WHbbLVl+551J+/ffT8q3356U64wcGXHUUR+Vf/ObiH79PiqPGBFx0kkflX/2s4ivfvWj8vXXR3z968nznj3zx7LVVh+1Hzo04vzzPypfdlnERRd9VP7Od5JHnYsuStrUOf/8pI86Z58d8cMfflQeODDimmsKx/KxjyXbUOekk5JtrNOvXzIGdY46KhmjOkcemYxhRDKmTc29xubKhsy9I49M5lBEMqeOPDKZYxFNz72dd84fy8c/nixv6dz7y1+S8ltvJeUxY5Lyu+8m5cbmXqH3p1Onj8a6JXMvInnfBw78qPzDHybzo06huVcolq233rC5V+frX09irPPVrzY997p2LTxfWjr3iv3ca2zutnTuFfu5161b459zLZl7EcV97hWaL507f9S+uXOvzoZ87nXpUvj92ZC5V8znXmNzpRR/c1vyubfLLo1/zhX7N7cln3s5c6VUyVwpDrPuCuTeO2tBWpe3TUSsBWqBHRpZtzl9AiBpsKQqSet+J8n8lt3OqyRWF7i1yOulv3VHk5YWuLXIgubf5qZkCm3/u82/nVfJFIrlvebfzqsk2tJcebPALZQWtew2NyVRaPtXtex2XiVRKJYVLbudV0ksz387r1aZL21p7i7JfzuvNvU5t7K09+BsUm3+23lt9nPljfy382pTn3PFKDYbBE4Dbs0pnwX8ukGbWUCPnPJcYEfg18CZOfV/SPtrss98j74NM/9NrdB/Zo7FsbTVOByLY3Es5RVLW4nDsTQrlra0Z24hsFtOuUdal7eNpPZAF2BpI+s2p8/COneG4S27RUdJDB+evLZjcSxZicOxOBbHUl6xtJU4HEvLYilWsdkgyTlw80guYKi7WKFPgzYXs+4FEPekz/uw7gUQ80gufmiyz3yPvnVZdmuc0FjnrruSGCTH4liyEYdjcSyOpbxiaStxOJYmYynVnrmS3M5L0gnAiDQRuy0ihku6BqiKiPGSOgF3Ap8GlgEDImJeuu4w4JvAWmBIRNxfqM+m4vDtvMzMzCwrSnU7L9+b1czMzKwV+N6sZmZmZuZkzszMzCzLnMyZmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLMCdzZmZmZhnmZM7MzMwsw5zMmZmZmWWYkzkzMzOzDHMyZ2ZmZpZhTubMzMzMMszJnJmZmVmGOZkzMzMzyzAnc2ZmZmYZ5mTOzMzMLMOKSuYkbS9pkqSX059dC7QblLZ5WdKgtK6zpH9IekHSbEk35LQ/W1KNpOnp47xi4jQzMzMrV8XumRsKTImIvYApaXkdkrYHrgI+CxwCXJWT9P08IvYFPg18XtLxOauOjYiD0setRcZpZmZmVpaKTeb6A3ekz+8ATsnT5jhgUkQsi4jlwCSgX0SsjIiHASLifeAZoEeR8ZiZmZltVopN5naOiEXp88XAznna7ArMzykvSOvqSdoOOIlk716dUyXNkDRO0m6FApA0WFKVpKqampoN2ggzMzOzrGoymZM0WdKsPI/+ue0iIoBoaQCS2gN3AzdHxLy0egLQKyIOINmTd0eh9SNiZERURERFt27dWvryZmZmZpnWvqkGEXF0oWWS3pDUPSIWSeoOvJmn2UKgMqfcA5iaUx4JvBwRI3Jec2nO8luBnzYVp5mZmdnmqNjDrOOBQenzQcB9edpMBI6V1DW98OHYtA5J1wFdgCG5K6SJYZ2TgeeLjNPMzMysLBWbzN0AHCPpZeDotIykCkm3AkTEMuBaYFr6uCYilknqAQwDegPPNPgKkkvTryt5FrgUOLvIOM3MzMzKkpJT3cpDRUVFVFVVtXYYZmZmZk2SVB0RFcX24ztAmJmZmWWYkzkzMzOzDHMyZ2ZmZpZhTubMzMzMMszJnJmZmVmGOZkzMzMzyzAnc2ZmZmYZ5mTOzMzMLMOczJmZmZllmJM5MzMzswxzMmdmZmaWYU7mzMzMzDLMyZyZmZlZhjmZMzMzM8swJ3NmZmZmGeZkzszMzCzDnMyZmZmZZZiTOTMzM7MMKyqZk7S9pEmSXk5/di3QblDa5mVJg3Lqp0p6UdL09LFTWt9R0lhJcyQ9JalXMXGamZmZlati98wNBaZExF7AlLS8DknbA1cBnwUOAa5qkPQNjIiD0sebad25wPKI+CRwE3BjkXGamZmZlaVik7n+wB3p8zuAU/K0OQ6YFBHLImI5MAno14J+xwFHSVKRsZqZmZmVnWKTuZ0jYlH6fDGwc542uwLzc8oL0ro6t6eHWH+Yk7DVrxMRa4FaYId8AUgaLKlKUlVNTU0Rm2JmZmaWPe2baiBpMrBLnkXDcgsREZKiha8/MCIWStoG+AtwFvDHlnQQESOBkQAVFRUtfX0zMzOzTGsymYuIowstk/SGpO4RsUhSd+DNPM0WApU55R7A1LTvhenPdyT9ieScuj+m6+wGLJDUHugCLG3OBpmZmZltToo9zDoeqLs6dRBwX542E4FjJXVNL3w4Fpgoqb2kHQEkdQBOBGbl6fc04KGI8F43MzMzswaa3DPXhBuAeySdC7wGfA1AUgVwQUScFxHLJF0LTEvXuSat24okqesAtAMmA79P2/wBuFPSHGAZMKDIOM3MzMzKkspph1dFRUVUVVW1dhhmZmZmTZJUHREVxfbjO0CYmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLMCdzZmZmZhnmZM7MzMwsw5zMmZmZmWWYkzkzMzOzDHMyZ2ZmZpZhTubMzMzMMszJnJmZmVmGOZkzMzMzyzAnc2ZmZmYZ5mTOzMzMLMOczJmZmZllmJM5MzMzswwrKpmTtL2kSZJeTn92LdBuUNrmZUmD0rptJE3PeSyRNCJddrakmpxl5xUTp5mZmVm5KnbP3FBgSkTsBUxJy+uQtD1wFfBZ4BDgKkldI+KdiDio7gG8Bvw1Z9WxOctvLTJOMzMzs7JUbDLXH7gjfX4HcEqeNscBkyJiWUQsByYB/XIbSNob2An4d5HxmJmZmW1Wik3mdo6IRenzxcDOedrsCszPKS9I63ININkTFzl1p0qaIWmcpN2KjNPMzMysLLVvqoGkycAueRYNyy1EREiKPO2aYwBwVk55AnB3RKyW9C2SvX5fKhDfYGAwwO67776BL29mZmaWTU0mcxFxdKFlkt6Q1D0iFknqDryZp9lCoDKn3AOYmtPHgUD7iKjOec2lOe1vBX7aSHwjgZEAFRUVG5pMmpmZmWVSsYdZxwOD0ueDgPvytJkIHCupa3q167FpXZ0zgLtzV0gTwzonA88XGaeZmZlZWWpyz1wTbgDukXQuydWoXwOQVAFcEBHnRcQySdcC09J1romIZTl9fA04oUG/l0o6GVgLLAPOLjJOMzMzs7Kkda85yLaKioqoqqpq7TDMzMzMmiSpOiIqiu3Hd4AwMzMzyzAnc2ZmZmYZ5mTOzMzMLMOczJmZmZllmJM5MzMzswxzMmdmZmaWYU7mzMzMzDLMyZyZmZlZhjmZMzMzM8swJ3NmZmZmGeZkzszMzCzDnMyZmZmZZZiTOTMzM7MMczJnZmZmlmFO5szMzMwyzMmcmZmZWYY5mTMzMzPLMCdzZmZmZhnmZM7MzMwsw4pK5iRtL2mSpJfTn10LtHtA0luS/t6gfg9JT0maI2mspC3T+o5peU66vFcxcZqZmZmVq2L3zA0FpkTEXsCUtJzPz4Cz8tTfCNwUEZ8ElgPnpvXnAsvT+pvSdmZmZmbWQLHJXH/gjvT5HcAp+RpFxBTgndw6SQK+BIzLs35uv+OAo9L2ZmZmZpajfZHr7xwRi9Lni4GdW7DuDsBbEbE2LS8Adk2f7wrMB4iItZJq0/ZLGnYiaTAwOC2uljSrZZuwWdiRPGNnHpc8PCb5eVzy87jk53FZn8ckv31K0UmTyZykycAueRYNyy1EREiKUgTVEhExEhgJIKkqIio2dQxtncclP4/L+jwm+Xlc8vO45OdxWZ/HJD9JVaXop8lkLiKObiSINyR1j4hFkroDb7bgtZcC20lqn+6d6wEsTJctBHYDFkhqD3RJ25uZmZlZjmLPmRsPDEqfDwLua+6KERHAw8BpedbP7fc04KG0vZmZmZnlKDaZuwE4RtLLwNFpGUkVkm6tayTp38CfSS5kWCDpuHTR/wW+K2kOyTlxf0jr/wDskNZ/l8JXyTY0ssjtKVcel/w8LuvzmOTnccnP45Kfx2V9HpP8SjIu8g4vMzMzs+zyHSDMzMzMMszJnJmZmVmGZSaZk9RP0ovpLb7WO4eusVuASboyrX8x53y9zGvGmHxX0nOSZkiaIqlnzrIPJE1PH+M3beQbVzPG5WxJNTnbf17OskHp7eleljSo4bpZ1oxxuSlnTF6S9FbOsrKcL5Juk/Rmoe+nVOLmdMxmSDo4Z1k5z5WmxmVgOh4zJT0u6cCcZa+m9dNL9bULbUUzxqVSUm3O78qPcpY1+vuXVc0YkytyxmNW+lmyfbqsnOfKbpIeTv8Gz5b0nTxtSvf5EhFt/gG0A+YCewJbAs8CvRu0uQi4JX0+ABibPu+dtu8I7JH20661t2kTjckXgc7p8wvrxiQtr2jtbWjFcTkb+HWedbcH5qU/u6bPu7b2Nm2qcWnQ/hLgts1gvhwBHAzMKrD8BOB+QMChwFPlPleaOS6fq9te4Pi6cUnLrwI7tvY2tNK4VAJ/z1Pfot+/LD2aGpMGbU8i+XaKzWGudAcOTp9vA7yU529RyT5fsrJn7hBgTkTMi4j3gTEkt/zKVegWYP2BMRGxOiJeAeak/WVdk2MSEQ9HxMq0+CTJd/mVu+bMlUKOAyZFxLKIWA5MAvptpDg3tZaOyxnA3ZskslYUEY8Ayxpp0h/4YySeJPluzO6U91xpclwi4vF0u2Hz+WxpznwppJjPpTathWOyWXyuAETEooh4Jn3+DvA8H93lqk7JPl+ykszV394rlXvrr/XaRPIlxHW3AGvOulnU0u06l+Q/gDqdJFVJelJS3nvqZlRzx+XUdLf2OEm7tXDdLGr2tqWH4/cAHsqpLtf50pRC41bOc6WlGn62BPCgpGolt1vc3Bwm6VlJ90vqk9Zt9vNFUmeShOQvOdWbxVxRctrXp4GnGiwq2edLsfdmtQyQdCZQARyZU90zIhZK2hN4SNLMiJjbOhFuchOAuyNitaRvkezR/VIrx9SWDADGRcQHOXWb83yxAiR9kSSZ+0JO9RfSubITMEnSC+nem83BMyS/KysknQD8DdirlWNqK04CHouI3L14ZT9XJG1NksAOiYi3N9brZGXPXN3tverk3vprvTZa9xZgzVk3i5q1XZKOJrmP7skRsbquPiIWpj/nAVNJ/msoB02OS0QszRmLW4G+zV03w1qybQNocCikjOdLUwqNWznPlWaRdADJ70//iKi/3WLOXHkTuJfyOK2lWSLi7YhYkT7/J9BB0o54vkDjnytlOVckdSBJ5EZHxF/zNCnd50trnyTYzBMJ25OcALgHH5082qdBm4tZ9wKIe9LnfVj3Aoh5lMcFEM0Zk0+TnHS7V4P6rkDH9PmOwMuUz8m4zRmX7jnPvwI8mT7fHnglHZ+u6fPtW3ubNtW4pO32JTkpWZvDfEm3qReFT2j/MuueoPx0uc+VZo7L7iTnH3+uQf1WwDY5zx8H+rX2tmzCcdml7neHJDF5PZ07zfr9y+qjsTFJl3chOa9uq81lrqTv+x+BEY20KdnnSyYOs0bEWknfBiaSXBV0W0TMlnQNUBUR40luAXankluALSNJ6Ejb3QM8B6wFLo51Dx9lUjPH5GfA1sCfk2tBeD0iTgY+BfxO0ocke2dviIjnWmVDSqyZ43KppJNJ5sMykqtbiYhlkq4FpqXdXRPrHhLIrGaOCyS/N2Mi/URJle18kXQ3yRWIO0paAFwFdACIiFuAf5JccTYHWAmcky4r27kCzRqXH5Gck/zb9LNlbURUADsD96Z17YE/RcQDm3wDNpJmjMtpwIWS1gLvAQPS36W8v3+tsAkl14wxgeSf5gcj4t2cVct6rgCfB84CZkqantZ9n+QfoZJ/vvh2XmZmZmYZlpVz5szMzMwsDydzZmZmZhnmZM7MzMwsw5zMmZmZmWWYkzkzMzOzDHMyZ2ZmZpZhTubMzMzMMuz/A3Tbz2RUY1llAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "HLYv61zCP7Cx" }, "source": [ "## 2-step Abysmal Butler Method \n", "The 2-step Abysmal Butler difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(4(t_i-w_i)-3(t_{i-1}-w_{i-1})), \\end{equation}\n", "for $i=1,...N.$\n", "For $i=0$ the system of difference equation is:\n", "\n", "\\begin{equation}w_{1} = w_{0} + \\frac{h}{2}(4(t_0-w_0)-3(t_{-1}-w_{-1})) \\end{equation}\n", "this is not solvable as $w_{-1}$ is unknown.\n", "\n", "For $i=1$ the difference equation is:\n", "\\begin{equation}w_{2} = w_{1} + \\frac{h}{2}(4(t_1-w_1)-3(t_{0}-w_{0})) \\end{equation}\n", "this is not solvable a $w_{1}$ is unknown. $w_1$ can be approximated using a one step method. Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.\\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "opfrOTvUP7Cx" }, "source": [ "### Initial conditions\n", "IC=1\n", "w=np.zeros(len(t))\n", "w[0]=IC\n", "\n", "w2=np.zeros(len(t2))\n", "w2[0]=IC\n", "w2[1]=(IC+1)*np.exp(-t2[1])+t2[1]-1\n", "\n" ], "execution_count": 5, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "cwkzmvkmP7C0" }, "source": [ "### Loop\n" ] }, { "cell_type": "code", "metadata": { "id": "wNl1YjAIP7C0" }, "source": [ "\n", "# Fine Mesh\n", "for k in range (1,N2):\n", " w2[k+1]=(w2[k]+h2/2.0*(4*myfun_ty(t2[k],w2[k])-3*myfun_ty(t2[k-1],w2[k-1]))) \n", " \n", " \n", "w[1]=w2[2]\n", "\n", "# Coarse Mesh\n", "for k in range (1,N):\n", " w[k+1]=(w[k]+h/2.0*(4*myfun_ty(t[k],w[k])-3*myfun_ty(t[k-1],w[k-1]))) \n" ], "execution_count": 6, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Pxr1ZksCP7C2" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "05okTgOJP7C3" }, "source": [ "def plotting(t,w,t2,w2):\n", " y=(2)*np.exp(-t2)+t2-1\n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,w,'^:',color='red',label='Coarse Mesh (N)')\n", " plt.plot(t2,w2, 'v-',color='green',label='Fine Mesh (2N)')\n", "\n", "\n", " plt.plot(t2,y, 'o-',color='black',label='Exact?')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.title('Abysmal Butler')\n", " plt.show " ], "execution_count": 7, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "1XA5KhK-P7C5" }, "source": [ "The plot below shows the Abysmal-Butler approximation for a low N (red) and $N\\times2$ (green) and the exact solution (black) of the intial value problem" ] }, { "cell_type": "code", "metadata": { "id": "yn18L9k2P7C5", "outputId": "291f63ef-e8e8-409b-f5af-b9d8b3cd2185", "colab": { "base_uri": "https://localhost:8080/", "height": 295 } }, "source": [ "plotting(t,w,t2,w2)" ], "execution_count": 8, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEWCAYAAABCPBKqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3zOV/vA8c/JnSV2i1JpEtReQSiqKLVVUTUataqhdPDr01YbikZadFCt2qvEVqO0RZXaJRF7zxhROxKReZ/fH99MGRIZd8b1fl73K/f5zutO88iVc873OkprjRBCCCGEeDJWlg5ACCGEECI3k2RKCCGEECIDJJkSQgghhMgASaaEEEIIITJAkikhhBBCiAyQZEoIIYQQIgMkmRJCPDGl1Hyl1DhLx5ESpVRzpdSVbLxfP6XUzuy6nxAiZ5BkSgjxWEqpbUqpu0opO0vHkpmUUlop9UApFaKUuqWUWqKUKpbGc11izrfO6jiFEDmbJFNCiFQppVyAlwANdLJoMFmjtta6EFAeKA6MsUQQSimTJe4rhMg4SaaEEI/TB9gLzAf6JrO/hFJqs1IqWCn1j1LKGUApNVUp9V3CA5VS65RSw2Pef6qUuhpz3imlVMuY7WOUUiuUUoti9h1RSlVSSn2mlLqhlLqslGqd4Jr9lVInYo49r5Qa9CQfUmt9H1gHVEtw7YtKqVcStMcopRbFNLfHfL0X07PV6NFrKqWqxHxv7sR8xu4J9s1XSk1TSv2ulHoAvPwkcQshLE+SKSHE4/QBfGJebZRSzzyy3x3wAkoAB2OOA1gA9FJKWQEopUoArwCLlVKVgfeA+lrrwkAb4GKCa74KLMToKfIHNmL8e1UW+BKYkeDYG0BHoAjQH5iklKqb3g+plCoOdMZIHNOiaczXYlrrQlrrPY9cryCwGVgMlAJ6Aj8rpaolOOxNwBsoDMhcKyFyKUmmhBApUko1AZyB5VprP+AcRgKQ0Aat9XatdTjgCTRSSj2ntd4HBAEtY47rCWzTWv8HRAN2QDWllI3W+qLW+lyCa+7QWm/UWkcBK4CSwHitdSSwFHCJnduktd6gtT6nDf8AmzCGJdPqgFLqHnALcCJxopYRHYGLWut5WusorbU/sAp4I8Exa7XWu7TWZq11WCbdVwiRzSSZEkKkpi+wSWt9K6a9mKRDfZdj32itQ4A7wLMxmxYAvWPe98bobUJrfRYYhjE/6YZSaqlSKvYcgP8SvH8I3NJaRydoAxQCUEq1U0rtjRlKuwe0x+glS6u6WutigD0wDdihlLJPx/kpcQZeUErdi31h9OKVTnDM5eRPFULkJpJMCSGSpZQqAHQHmimlriulrgPDgdpKqdoJDn0uwTmFgKeAazGbFgGvxRxfFVgTe6zWerHWOrbnSwMTniBGO4zenm+BZ2KSot8Bld5rxfR6zQbKATViNj8AHBIcljAR0o+55GXgH611sQSvQlrrd9NxDSFELiDJlBAiJZ0xhuOqAa4xr6rADox5VLHaK6WaKKVsMeZO7dVaXwbQWl8B9mP0SK3SWj8EUEpVVkq1iEmGwjB6m8xPEKMtxnDhTSBKKdUOaJ36KcmLeZquf0ws52M2HwR6KqVslFJuQLcEp9yMibl8CpdcD1RSSr0Vc76NUqq+Uqrqk8QnhMi5JJkSQqSkLzBPax2gtb4e+wJ+AtwT1FdaDIzGGN6rR/ywXqwFQE1ihvhi2AHjMeYpXceYoP1ZegPUWgcDHwDLgbsY87nWpfMyh5RSITHn9wW6aK3vxOwbBVSI2TcW47PG3jsUY/L4rphhvIbJxNYaY67YNYzPOQHjswsh8hCltfQyCyGyjlKqKcZwn7OWf3CEEHmQ9EwJIbKMUsoG+BCYLYmUECKvkmRKCJElYuYG3QPKAJMtHI4QQmQZGeYTQgghhMgA6ZkSQgghhMgAi612XqJECe3i4mKp2wshhBBCpJmfn98trXXJ5PZZLJlycXHB19fXUrcXQgghhEgzpdSllPbJMJ8QQgghRAZIMiWEEEIIkQGSTAkhhBBCZIDF5kwlJzIykitXrhAWFmbpUEQmsre3x9HRERsbG0uHIoQQQmS6HJVMXblyhcKFC+Pi4oJS6V70XeRAWmtu377NlStXKFeunKXDEUIIITJdjhrmCwsL4+mnn5ZEKg9RSvH0009Lb6MQQohM5+Pjg4uLC1ZWVri4uODj42OROHJUzxQgiVQeJP9NhRBCZDYfHx88PDwIDQ0F4NKlS3h4eADg7u6erbHkqJ4pIYQQQoi08PT0jEukYoWGhuLp6ZntsUgyJYQQQohcJyAgIF3bs1LuT6YCA6FZM7h+PVMud/36dXr27EmFChWoV68e7du35/Tp05ly7Sc1f/58lFL89ddfcdvWrFmDUoqVK1em+3pjxozh22+/fexxWmtatGjB/fv3AWO47qOPPorb/+233zJmzBgAfvrpJ+bOnZvuWIQQQoj00lpTuHDhZPc5OTllczR5IZny8oKdO42vGaS1pkuXLjRv3pxz587h5+fH119/zX///Zfh65rN5gxdo2bNmixdujSuvWTJEmrXrp2haz7O77//Tu3atSlSpAgAdnZ2/Prrr9y6dSvJsQMGDODHH3/M0niEEEIIrTXDhw/n/v37WFsnnvrt4OCAt7d3tseUs5Op5s1h/nzjfWSk0V60yGiHhkKjRjB7NpjNMHcuNG4Mv/5q7L91yzj+t9+Mdhp6rrZu3YqNjQ2DBw+O21a7dm1eeukltNZ8/PHH1KhRg5o1a7Js2TIAQkJCaNmyJXXr1qVmzZqsXbsWgIsXL1K5cmX69OlDjRo1uHz5Mv369Ys7f9KkSQCcO3eOtm3bUq9ePV566SVOnjyZbGwvvfQS+/btIzIykpCQEM6ePYurq2vcfj8/P5o1a0a9evVo06YNgYGBAEyZMoVq1apRq1YtevbsGXf88ePHad68OeXLl2fKlCnJ3tPHx4fXXnstrm1tbY2Hh0dc7Ak5ODjg4uLCvn37Hvt9FkIIIZ6E2Wxm6NCh/PDDD3z44YfMnz8f50KFUIBzoULMnDkz2yefQw58mi9dLl0CrY33ZrPRzoCjR49Sr169ZPf9+uuvHDx4kEOHDnHr1i3q169P06ZNKVmyJKtXr6ZIkSLcunWLhg0b0qlTJwDOnDnDggULaNiwIX5+fly9epWjR48CcO/ePQA8PDyYPn06FStW5N9//2XIkCH8/fffSe6vlOKVV15h48aNBAUF0alTJy5cuAAYxU7ff/991q5dS8mSJVm2bBmenp7MnTuX8ePHc+HCBezs7OLuCXDy5Em2bt1KcHAwlStX5t13301SVHPXrl3MmDEj0bahQ4dSq1YtPvnkkyQxurm5sWPHDho0aJDWb7kQQgiRJtHR0QwaNIg5c+bwySefMH78eNSFC7hHRcUeAC1bWiS2nJ1MbdsW/97GJnE7KAju3oXYb2JEhNFu3NholyiR+PjSpTMUys6dO+nVqxcmk4lnnnmGZs2asX//ftq1a8fnn3/O9u3bsbKy4urVq3HDgs7OzjRs2BCA8uXLc/78ed5//306dOhA69atCQkJYffu3bzxxhtx9wkPD08xhp49ezJlyhSCgoL47rvv+OqrrwA4deoUR48epVWrVoDxA1emTBkAatWqhbu7O507d6Zz585x1+rQoQN2dnbY2dlRqlQp/vvvPxwdHRPd786dO0nGpIsUKUKfPn2YMmUKBQoUSLSvVKlSKfasCSGEEE8qKiqKAQMGsHDhQr744gvGjBmDWrsWevUyOlPASKa8vGDq1GyPL2cP86XGyyv+Gxgr9hv5hKpXr46fn1+6zvHx8eHmzZv4+flx8OBBnnnmmbgClQULFow7rnjx4hw6dIjmzZszffp0Bg4ciNlsplixYhw8eDDudeLEiRTv1aBBA44cOcKtW7eoVKlS3HatNdWrV4+7xpEjR9i0aRMAGzZsYOjQoRw4cID69esTFZN82tnZxZ1vMpnitidkbW2d7FyvYcOGMWfOHB48eJBoe1hYWJIESwghhMiIyMhI3N3dWbhwIePGjWPsqFFG/cKyZY2OlIgI48CICJg3L9MeSEuP3JtM7dkT/w2MFREBu3c/8SVbtGhBeHg4M2fOjNt2+PBhduzYwUsvvcSyZcuIjo7m5s2bbN++nQYNGhAUFESpUqWwsbFh69atXEphqPHWrVuYzWZef/11xo0bx4EDByhSpAjlypVjxYoVgJEUHTp0KNUYx48fH9cjFaty5crcvHmTPXv2AMYP3rFjxzCbzVy+fJmXX36ZCRMmEBQUREhISJq/H5UrV+b8+fNJtj/11FN0796dOXPmJNp++vRpatSokebrCyGEEKmJiIigR48eLF++nG8mTsTz9GkYMMDYOW8ePDIBPaOdKk8q9yZT/v7GfKlHX/7+T3xJpRSrV6/mr7/+okKFClSvXp3PPvuM0qVL06VLF2rVqkXt2rVp0aIFEydOpHTp0ri7u+Pr60vNmjX55ZdfqFKlSrLXvnr1Ks2bN8fV1ZXevXvz9ddfA0bP1pw5c6hduzbVq1ePm8Ceknbt2vHyyy8n2mZra8vKlSv59NNPqV27Nq6uruzevZvo6Gh69+5NzZo1qVOnDh988AHFihVL8/ejQ4cObEs4VJrARx99lOSpvl27dsUNNQohhBAZERYWRteuXVm9ejVTpkzhfx9/DBUrQoUKxu/7LOhUeVJKx07gzmZubm7a19c30bYTJ05QtWpVi8QjkgoMDKRPnz5s3rz5scf6+/vz/fffs3DhwmT3y39bIYQQaRUaGkqXLl3YtGkTM8qWxWP1aqhf36IxKaX8tNZuye3LvT1TIsuVKVOGd955J65oZ2pu3bqFlwW6VoUQQuQtISEhdOjQgc2bNzN36lQ8nJ3hkTm6OU3OfppPWFz37t3TdJwM7wkhhMio+/fv075ePfacPcvChQtx790bhgyxdFiPJcmUEEIIISzu3r17tG3bFr8LF1jq6sobr75q6ZDSTIb5hBBCCGERPj4+uDg5YaUUpUqUYP/+/axcuZI3DhyAokUtHV6aSc+UEEIIIbKdj48PHh4ehIaGAhAZHY2dnR0hDx6AUhaOLn2kZ0oIIYQQ2cvXF8933olLpGKFh4fj6elpoaCeXK5NpurMqIMaq5K86syok6HrmkwmXF1d414XL16kcewSNRk0ZswYlFKcPXs2btvkyZNRSvFomYi06NevHytXrnzscQ8fPqRZs2ZER0dz8OBBGjVqRPXq1alVq1bcgs0AzZs3x80t/qlPX19fmjdvDsCRI0fo169fumMUQgghkggKIuDhw2R3BQQEZHMwGZdrk6lGjo2wNdkm2mZrsqWxY8YSnwIFCiRa3sXFxYXdmVgArGbNmixdujSuvWLFCqpXr55p10/O3Llz6dq1KyaTCQcHB3755ReOHTvGn3/+ybBhwxItgHzjxg3++OOPZOO+cuVKrvwhF0IIYWFaw7Rpxgs48eyzWJlMyR7q5OSUnZFlihw7Z2rYn8M4eP1givvDo8KJMideTy7KHIX/dX+az2+e7DmupV2Z3HZyumMpVKgQISEhbNu2jTFjxlCiRAmOHj1KvXr1WLRoEUop/Pz8+L//+z9CQkIoUaIE8+fPj1tsOKHOnTuzdu1aRo4cyblz5yhatCg2NjZx+zdt2sTo0aMJDw+nQoUKzJs3j0KFCjFixAjWrVuHtbU1rVu35ttvvwVg+/btfP/991y/fp2JEyfSrVu3JPf08fFh8eLFAInW9Hv22WcpVaoUN2/ejKuM/vHHH+Pt7U27du2SXOfVV19l6dKlfPLJJ+n+HgohhMjHlII//gCl8H/hBVq3aUPBggWJiIiIW88WwMHBAW9vbwsG+mRybc+UnbUdzxR8BoUxSU2hKF2wdJLeqvR6+PBh3BBfly5dkuz39/dn8uTJHD9+nPPnz7Nr1y4iIyN5//33WblyJX5+fgwYMCDFMd8iRYrw3HPPcfToUZYuXUqPHj3i9t26dYtx48bx119/ceDAAdzc3Pj++++5ffs2q1ev5tixYxw+fJiRI0fGnRMYGMjOnTtZv349I0aMSHK/iIgIzp8/j4uLS5J9+/btIyIiggoVKsRta9SoEba2tmzdujXJ8W5ubuzYsSPV758QQggBwO3b8P77cOOG0V68mD2ffsrLLVpQoEABfH19mT17Ns7OziilcHZ2ZubMmbi7u1s27ieQY3um0tKDFBgcSPkp5QmLCsPe2h6/QX6ULlQ6Q/eNHeZLSYMGDXB0dASIm1NVrFgxjh49Gle4Mjo6OtleqVg9e/Zk6dKlbNy4kS1btjBv3jwA9u7dy/Hjx3nxxRcBIxFq1KgRRYsWxd7enrfffpuOHTvSsWPHuGt17twZKysrqlWrxn///ZfkXrdu3Up2Pb7AwEDeeustFixYgJVV4px65MiRjBs3jgkTJiTaXqpUKa5du5bi5xJCCCHi3LwJc+dCkybQowdb/v2X1157jTJlyrBlyxacnJyoWLFirkyeHvXYniml1Fyl1A2l1NEU9rsrpQ4rpY4opXYrpWpnfpjJK1O4DP1d+2OlrOjv2j/DiVRa2NnZxb03mUxERUWhtaZ69epx86yOHDnCpk2bUrxGx44dWbhwIU5OThQpUiRuu9aaVq1axV3n+PHjzJkzB2tra/bt20e3bt1Yv349bdu2TTae5NZZLFCgQKIuVDAqzHbo0AFvb28aNmyY5JwWLVrw8OFD9u7dm2h7WFgYBQoUSOW7I4QQIt8JDIRmzeD6ddi/H77/3thepQoEBECPHvz222906NCBcuXKsX379lw5Lyo1aRnmmw+0TWX/BaCZ1rom4AXMzIS40mxU01E0cWrCqGajsvO2iVSuXJmbN2+yZ88eACIjIzl27FiKxzs4ODBhwoQkQ4ENGzZk165dcU/7PXjwgNOnTxMSEkJQUBDt27dn0qRJHDp0KM2xFS9enOjo6LiEKiIigi5dutCnT59k51fFGjlyJBMnTky07fTp09SoUSPN9xZCCJEPeHnBzp3GVx8fI5kKDjb2Pf00S5cupWvXrtSsWZNt27alOnKTWz12mE9rvV0p5ZLK/oSPuu0FHDMeVtqVKVyGf/r9k523TMLW1paVK1fywQcfEBQURFRUFMOGDUv1Kb2ePXsm2VayZEnmz59Pr169CA8PB2DcuHEULlyY1157jbCwMLTWfB+b9adR69at2blzJ6+88grLly9n+/bt3L59m/nz5wMwf/58XF1dE53Tvn17SpYsmWjb1q1b6dChQ7ruLYQQIg+7cgVmzwazGebNg8OH4csvoXBhAGbPno2HhwdNmjRh/fr1iUZj8hKV3NBQkoOMZGq91jrVbgml1P+AKlrrgSns9wA8AJycnOpdunQp0f4TJ05QtWrVNAUu0u7AgQNMmjSJhQsXPvE1wsPDadasGTt37sTaOv1T7eS/rRBC5EFvv23MiwKwtYWBA2HqVMCoozh8+HDatGnDr7/+ioODgwUDzTillJ/W2i25fZn2NJ9S6mXgbeDTlI7RWs/UWrtprd0e7fUQWadu3bq8/PLLREdHP/E1AgICGD9+/BMlUkIIIfKQ48dh1Ci4dg1iyu4AEBEB8+ahAwPx8vJi+PDhdO3albVr1+b6ROpxMuU3o1KqFjAbaKe1vp0Z1xSZa8CAARk6v2LFilSsWDGTohFCCJFrbdkCP/4IFy8aw3sJ6KgoPm3fnm8OHuStt95i7ty5+eKP8Az3TCmlnIBfgbe01qczHpIQQgghcoyICJgwAf7802gPHgxnz8LRo8a+GGZgSGQk3xw8yLvvvsv8+fPzRSIFaSuNsATYA1RWSl1RSr2tlBqslBocc8gXwNPAz0qpg0qp9C8yJ4QQQoicSSljXtTGjUbbxgZKlAB/f3wWLcLF2RkrpShSsCDTgU8++YSpU6cmqWGYl6Xlab5ej9k/EEh2wrkQQgghciFfX5g0CebPN5KnvXuhePFEh/j4+ODh4UFoaChglPOxsbGhVq1aKKUsELTl5J+0MY1MJlPccjKurq6MHz8+06598OBBfv/997j29OnTqV69OpUqVWLMmDGZdh8hhBAiQwIDYetWYzgPkiRSAJ6ennGJVKzIyMgUl1PLbHVm1EGNVUledWbUyZb7J5SrkykfHx9cXFywsrLCxcUFHx+fDF8zdjmZ2Fdy6909qUeTqeeffx5/f3+OHDnCggULuHLlSqbdSwghhEizyEj49FOYMcNod+xoJFKplLQJCAhI1/bM1sixUZL1eG1NtjR2bJwt908o1yZTsd2Lly5dQmvNpUuX8PDwyJSE6lFBQUFUrlyZU6dOAdCrVy9mzZoFwLvvvoubmxvVq1dn9OjRcefs37+fxo0bU7t2bRo0aEBQUBBffPEFy5Ytw9XVlWXLlvHKK69ga2uL1pqoqChsbTO2SLMQQgjxRKytwd8fYn7PoRSkUs7g6tWrmEymZPdl11IxI14cgVknfprQpEwWWRElx06zHzZsWKoLDu/duzeuSnis0NBQ3n777bhE51Gurq5Mnpz6AsoPHz5MVA38s88+o0ePHvz000/069ePDz/8kLt37/LOO+8A4O3tzVNPPUV0dDQtW7bk8OHDVKlShR49erBs2TLq16/P/fv3cXBw4Msvv8TX15effvop0T09PDzo2bMnpUqVSjU2IYQQItPs2QOffALr1hnDeBs2GPOjHuP06dO0bt0ak8mEyWRK9LvYwcEBb2/vrIwagO2XtjNkwxCizFEoFBqNrck229bpfVSOTaYe59FE6nHb0yp2mO9RrVq1YsWKFQwdOjTR2njLly9n5syZREVFERgYyPHjx1FKUaZMGerXrw+Qavn8devWERgYGLe0ixBCCJGltI7vebp1Cy5fNpKpNCRSBw4coG3btmit2bVrFydPnsTT05OAgACcnJzw9vbG3d09y0K/HnKdjzd/zKLDi3Au6szcTnMZ8vsQwqLCLNYrBTk4mXpcD5KLiwuPLkcD4OzszLZt2zI9HrPZzIkTJ3BwcODu3bs4Ojpy4cIFvv32W/bv30/x4sXp169f3ILCaXX48GFat26drx4hFUIIYQFaGzWiihaFiROhdm04dgzS+Ptn27ZtdOrUieLFi7Np0yYqV65MvXr1sjR5ihVljmLa/mmM3DqSh5EP+bzJ53g29cTBxoH91/Yzw2+GxXqlIBfPmfL29k5Snj4ruxcnTZpE1apVWbx4Mf379ycyMpL79+9TsGBBihYtyn///ccff/wBQOXKlQkMDGT//v0ABAcHExUVReHChQmOXUk7RufOnenUqVOWxCyEEEIQu5SYUmAyGV9jpTGRWrt2LW3btsXR0ZFdu3ZRuXLlLAg0eXsu76H+rPp88OcHvFD2BY4OOYp3S28cbIwcYFTTUTRxamKxXinIxcmUu7s7M2fOxNnZGaUUzs7OzJw5M8MZcuycqdjXiBEjOHXqFLNnz+a7777jpZdeomnTpowbN47atWtTp04dqlSpwptvvsmLL74IgK2tLcuWLeP999+ndu3atGrVirCwMF5++WWOHz8eNwEdYOfOnfz7778Z/n4IIYQQSezZA5UqwZkzRnvqVKOaeTrMnz+f119/ndq1a7Njxw4cHR2zINCkbj64ydtr36bx3MbcfHCT5d2Ws7H3Rio9XSnRcWUKl+Gffv9YrFcKQGmtLXJjNzc37eubuFj6iRMnqJrKY5gi95L/tkIIkY0iIsDWFq5fh169YPJkY1gvnb7//ns++ugjWrVqxa+//kqhQoWyINjEos3RzD4wm8+2fEZwRDDDXhjGF82+oLBd4Sy/d2qUUn5aa7fk9uXYOVNCCCGESIPAQOjZE5Ytg9KloV8/uH8ffv3VaG/dmu5Laq3x9PTk66+/5o033mDhwoXY2dllath1ZtTh4PWkD3w5WDsQGhVKM+dmTG0/leqlqmfqfbNCrh3mE0IIIQTg5QU7dxpfweiBqlsXzObUz0tBdHQ0gwYN4uuvv2bQoEEsWbIk0xMpSL7oJoBGs6jLIrb23ZorEinIgcmUpYYdRdaR/6ZCCJFFAgONRYjNZpgzxxjWGz4cRo5M8+TyhMLDw+nZsyezZs3i888/Z9q0aSkW58yoUU1HYaUSx2hSJg4OOoh7Lfdctb5fjkqm7O3tuX37tvzyzUO01ty+fRt7e3tLhyKEEHlHeDhcvGj0RsX+zjSb43un0ijhsmxOTk7UrVuXlStX8t133+Ht7Z2lCc3FexcpbBs/D8rGyoZB9QZRqUSlVM7KmXLUBPTIyEiuXLmS7lpNImezt7fH0dERmzQUhBNCCJEGLVoYBTfPnIGEvzMLFIDz5425Uo8Ruyzbo4sVe3h4MCN2jb4scC34GiP+GsHCwwsp5VCKu2F3iTRHUsC6AOc/PG/Rp/JSk2smoNvY2FCuXDlLhyGEEELkLGYzrF0LnToZtaI+/hh+/jl+Lb1Y0dFG79TUqY+9pKenZ5JECmDjxo2ZFXUi4VHhTNo7iXHbxxFpjuSzJp/x+Uuf88nmTyxedDOjctQwnxBCCCGS8ccf0LUrrFljtNu1gytXjBIICUVEwO7dabpkQEBAurY/Ka01606to/rP1flsy2e0LN+S40OO81XLryhkWyhHFN3MqBzVMyWEEEIIjHlQmzdDaCh07mwkT+vWQYcO8cf4+2foFqVKleK///5Lst3JySlD103oxM0TDNs4jE3nNlG1RFU29t5I6wqtEx0TW3QzN5NkSgghhMiJvvzS+Nq5s/Fk3quvZtqlV61axa1bt1BKJXroK7OWZQsKC2LsP2P5cd+PFLQpyKQ2kxhafyg2prw5d1aSKSGEECInOHrUmO80ezYULgxLlkCpUpl+m59++okPPviAhg0b0rdvX77++msCAgJwcnLC29s7zcuypVR087kizxEWFcat0FsMrDuQcS3GUapg5n+OnESSKSGEEMKStDYWHw4Nhb//NpKqRo3guecy+TbxVc07derEkiVLcHBwYNCgQU90vUaOjTh+8zgR0fHzthSKy/cv8+JzL/KH+x/Ue7ZeZoWfo8kEdCGEEMISoqOhd28YFTPxukEDCAgwEqlMFhkZSf/+/fn666/x8PBg1apVODg4ZOiayRXd1Gimtp/Kjv478k0iBemq70cAACAASURBVJJMCSGEENkrONj4ajIZdaESFjUuUCDTbxcSEkKnTp1YsGABY8eOZfr06VhbZ3xgqrBdYaqXjF/uxaRMDKw7kCH1h+Sq6uWZQYb5hBBCiOyyZAkMGgTHjhnDeLNmZentbty4QYcOHfD392fWrFkMHDgww9eMNkez4NACRv49ksCQQKyUFWZtxtZki9fL6avAnldIz5QQQgiRlW7fhmvXjPcvvghvvgnZsCLEuXPnaNy4MceOHWPNmjWZkkj9df4v6s6sy9vr3sa5mDO7B+xmUL1BWCmrXF10M6OkZ0oIIYTIKuHhUKMGtGwJixaBkxNMn57lt/X19aV9+/aYzWb+/vtvGjZsmKHrnbh5go83f8yGMxtwKebC0teX0r16d5RSuBRz4djNY7m66GZGSc+UEEIIkVGBgdCsGVy/Dg8ewMqVxnY7O/juO/j882wLZePGjTRv3pyCBQuya9euDCVSNx/cZOiGodScVpMdATuY+MpETgw9QY8aPeLmRcUW3cyvvVKQhmRKKTVXKXVDKXU0hf1VlFJ7lFLhSqn/ZX6IQgghRA7n5QU7dxpff/oJ3ngjft28N9+EatWyJYxffvmFjh07UrFiRXbv3k3lypWf6DphUWFM3DWR5398nhl+MxjsNpiz75/l4xc/xt7a/vEXyGfSMsw3H/gJ+CWF/XeAD4DOmRSTEEIIkXsEBBiFNs1mmDcPDh+GJk3gCROZ9PDx8cHT05OAgACKFi3KvXv3aNmyJb/++itFihRJ9dyUim46F3UG4FLQJTpW6sjEVyZStWTVLIk/r3hsz5TWejtGwpTS/hta6/1AZGYGJoQQQuRoscuweHtDZMyvwOhomDTJmGiexXx8fPDw8ODSpUtorbl37x4mk4m33nrrsYkUGEU3bU22ibYpFJeCLlHMvhh/vfUXv/X6TRKpNMjWOVNKKQ+llK9SyvfmzZvZeWshhBAi88ycCa1aGU/p/ZJg4CYiwuidun49y0Pw9PQkNDQ00bbo6GhGjx6dpvNTKro5qfUk/Dz8aFm+ZabFmtdlazKltZ6ptXbTWruVLFkyO28thBBCZMyVK/E9UHZ2UKiQUb3cbE58XHS0MXcqiwUEBKRr+6NMViaeL/58fDum6OawRsMwWZkyJcb8Qp7mE0IIIR7n0CEoXx6WLTPaffvCmjVw4IDRG5VQRATs3p2l4Zw4cQIrq+R/hTs5OaV6bnB4MGO2jaHClAocv3kckzISp/xcdDOjJJkSQgghknPgAPz+u/G+Zk0YOdKYWJ6Qv78xd+rRl79/loW1detWGjduTMGCBbG3T/xknYODA97e3smeFxEdwY///kiFKRUY+89Y2j7fluNDj+NRzyPfF93MKKVjJ9CldIBSS4DmQAngP2A0YAOgtZ6ulCoN+AJFADMQAlTTWt9P7bpubm7a19c3o/ELIYQQWaN5c7h1C44cgRyy1twvv/zCwIEDqVSpEhs2bGDnzp1xT/M5OTnh7e2Nu7t7onPM2szSo0sZ+fdILty7QHOX5kx4ZQINyjYAIDA4kJ6rerKs2zJJplKhlPLTWrslu+9xyVRWkWRKCCFEjrJvH3zxBSxfDkWKwNmzULIkFC1q6cjQWjN27FjGjh1LixYtWLVqFcWKFXvsOZvObWLElhEcvH6Q2s/UZvwr42lToU2+W4g4M6SWTMlyMkIIIfKvyEhjyZdChcDKyii0ee4c1KkDzz//+POzQUREBAMHDmThwoX069ePGTNmYGtrm+o5+6/uZ8SWEfx94W9cirmwqMsietXsleTpPZE5JJkSQgiRPz18aMyFev11mDAB3NyM3ihTznmS7e7du3Tt2pVt27bh5eWFp6dnXK9SSkU3i9kV4174PUo4lOCHtj8wqN4g7Kztsjv0fEWSKSGEEPlHUJCx7EuHDlCgAPTrZyRRsXJQInXhwgXat2/P+fPnWbRoUZK5UI0cG3H85nEiohM/TRgSEcIXTb/go8YfUcTu8cU7RcbJnCkhhBD5x/DhMHWqUWyzRAlLR5Oiffv28eqrrxIZGcnq1atp1qxZkmMCgwMpP6U8YVFhcdtMyoT/IH9qPlMzO8PNF1KbMyWDp0IIIfKuGzdgyBA4ftxof/QR7N2boxOp1atX07x5cwoWLMju3buTTaTuhd1juu90os3RcdtsrGwYVG+QJFIWIMmUEEKI3C8wEJo1i1/GJbZSuckES5fCv/8abUdHqFvXMjGmwMfHBxcXF6ysrHjqqafo2rUrtWrVYu/evVSpUiXRsSERIXy14yvK/VCOL7d/SavyrbAzGfOhrK2sGdVslCU+Qr4nc6aEEELkfl5exlwoLy/j6bzr12H9enj6abh8GQoWtHSEyYpdrDh2jb27d+9iMpkYNGgQpUqVijsuNDKUn/f/zIRdE7gVeouOlTryZfMvqVOmDkM2DGGG3wwpumlBMmdKCCFE7nb1qrHUS0SEManc09NYH2/UqBxTbDMlLi4uXLp0Kcl2Z2dnLl68SHhUODP9ZvLVzq+4HnKdVuVb4fWyFy84vhB3rBTdzB5SZ0oIIUTe1bdv/Pp40dHG5PKpUy0bUxqltljxTL+ZeG334sr9KzR1bsrybst5yfmlJMeWKVyGf/r9k9WhilTInCkhhBC5S0gIvPsurFplzJXatSt+X0QEzJsXP3cqB9uzZw+a5EeHdBHNoPWDcCziyOa3NrOt77ZkEymRM0gyJYQQIueLijIKagI4OBgJ1NmzxhwpsznxsdHRxvYcbOHChTRv3hz7IvZJx4hsoGDbgmx4cwO7B+zmlfKvyPIvOZwkU0IIIXK+3r2hVSsjUbKyAn9/+PRT2LMnfogvVkQE7N5tmTgfw2w289lnn9GnTx8aN27MPv99WL1mBbHL/xUF02smzsw6Q/uK7SWJyiVkzpQQQoic5/Rp+PZb+O47KFzYGNbr0SN+f2ylcn9/y8T3BEJCQujduzdr167l7YFv02RwE7r/3h1zTTPElIayNdkysM5AyhQuY9lgRbpIz5QQQoicISICgoON93fuwOLFcOCA0W7WDLp0yVHLvaRHQEAATZo04bfffqPn/3qytcZW+q/vj42VDdM7TMfe2h4wKphLrajcR5IpIYQQlvfggVHeYPx4o/3CC/GFOHO5PXv20KBBA06dPUWJgSVYWmgpRe2L8mv3Xzk4+CCD3AbR37U/VspKakXlUjLMJ4QQwjL27oVDh2DQIKOo5rvvQuPGxj6ljOG9XG7egnl4eHigC2ui+0VT27U2c5vOTTIfalTTURy7eUx6pXIpKdophBAi+0REgK2t8X7oUKO8waVLYGdn2bgyWXBYMF0GdWHLL1vAGRp+1JAv238pT+blYlK0UwghhOVt3QpvvAE7dkDVqjB2LEycmKsTqToz6nDw+sHEG8OBNcAJeLbZsyyYuYBXKr1iifBENpFkSgghRNbQGv7+G556CurUgRo1oEWL+P0lSlgutkzSyLERR/46QvTmaAgCCgMKCIYPvviAyWMmS09UPiAT0IUQQmRc7GTx69eNJAqMBYd79jTKGwCULAnLlxu9UnlE2QtliV4bk0gBBAP34d0P3uWHsT9IIpVPSM+UEEKIjPPygp074fXXjcWGN28Ge3vYtClPJU+xTt46yTe7vmHuyLkQmXT/72t+h8nZH5ewDEmmhBBCPLmQEGMtvHnzjGVd9u+HV1+FsDAjqapTx9IRZqo9l/cwYdcE1p5aa9SGup/8cSktYCzyJhnmE0IIkT5mszGEB7BtG3zwgbF2HhglDUqXNhKpPMKszaw/vZ6X5r1E47mN2X5pO6OajuJQn0PY29kne46Tk1M2RyksSZIpIYQQaXf3LlSoAD//bLRr1zZKHcQmUxERRi/V9euWizGTRERHsODgAmpNq8WrS17l0r1LTG4zmYDhAXR7uhttm7UlMjIS29hSDzEcHBzw9va2UNTCEmSYTwghROoWLoT79426UMWLQ8eO8fOgvv466fHR0cYcqqlTszfOJ5BsaQPg2cLPYqWsuHL/CjVL1WRhl4X0qN4DG5MNy5cvp3///hQrVoydO3dy7tw5PD09CQgIwMnJCW9vb9zd3S3waYSlPLZnSik1Vyl1Qyl1NIX9Sik1RSl1Vil1WClVN/PDTJ86M+qgXleoYgqlYr6+rqgzI2+N3QshRJaIjk68gPD69bB0aXz7xx+hbVvj/Z49Rm9UQhERsHt31seZCRo5NsLWZJtk+7Xga1QoXoENb27g0OBD9K7VGyus+PTTT+nRoweurq74+vrSsGFD3N3duXjxImazmYsXL0oilQ+lZZhvPtA2lf3tgIoxLw9gWsbDypinzzwNvxH/qGoQ8BuUOJP7a5oIIUSW++orqF8f/vvPaM+eDdu3J3+sv79RCuHRV8JkLAcb1XQUisTlC6yUFRt6bWBbv21xy77cuXOH9u3bM3HiRAYPHszWrVspU6aMhaIWOc1jkymt9XbgTiqHvAb8og17gWJKKYv+hJ1cfjLpo6qRMduFEEIkduKEsSZe7BJf7u6weDEUK2a0Cxc2JpbnIVprNp7dyIB1AwiPDo/bbmNlw+B6g2lfqX3ctiNHjlC/fn22bt3KzJkzmTZtWpJ5UiJ/y4wJ6GWBywnaV2K2JaGU8lBK+SqlfG/evJkJt07etSvXkt1+9crVLLunEELkGmYz/PmnMUQHUKaMMTR3757RLl8eunfP1cu8pCQ0MpQZvjOo/nN12vq05eD1g3zS+BPsTcZTedZW1okWG16+fDkNGzbk4cOH/PPPP7zzzjuWCl3kYNk6AV1rPROYCcZCx1l1H6dnn+XS1aSJk5N0yQoh8rOgICha1BiGGzgQXnwRGjUyeqDy+MLzV+5fYeq+qcw8MJM7D+9Qt0xdFnZZSPfq3bE12RIcEcwMvxn0d+1P6UKliY6OxtPTkwkTJtCoUSNWrVolw3oiRZmRTF0FnkvQdozZZjHeVargcfUqoQm2KQXed+/CjRtQqpTFYhNCCIsYNMhYYPjYMTCZjMrkFSpYOqos9++Vf5n872RWHFuBRtOlSheGNRzGi8+9mGipl1FNR3Hs5jFGNRvFnTt3ePPNN9m4cSMeHh5MmTIFuzzYSycyT2YkU+uA95RSS4EXgCCtdWAmXPeJud++DcCnCq5qsDFBZDRsKxuNex5YWFMIIZIVGGishbdsGVy9CpMmwaxZRgHN9u2hShWjHpSNDVSrZuloM0VKpQ2cizpTpnAZ9l7ZSxG7IgxrOIz3GryHSzGXJMf6+PjElTao+3ldoqKiuHfvHjNmzMDDwyMbPoXI7R6bTCmllgDNgRJKqSvAaMAGQGs9HfgdaA+cBUKB/lkVbJr5++MOxD6cGhYRRhnXMsw+c4/mm5fgXq8NLFkC772X5yZVCiHyKa3hww+N9fG8vKBrV6P36dQpcHWF116zdIRZopFjI47fPE5EdOLyDJeCLmFjsuHHdj/St3ZfCtsVTvZ8Hx8fPDw8CA01xjICA42+gNGjR0siJdJMaZ1lU5dS5ebmpn2zcYz+xMUT1KxTExSceG84FSf+CIcOQeXK2RaDEEJkKrMZHjwwnrbz8wM3N2N7gQJw9iyULGn0QuVhgcGBlPuhXKIn8qyUFfNfm497LXesVOrPWbm4uHDp0qUk252dnbl48WJmhytyMaWUn9baLbl9+WY5maouVZmzaA7R96NpuHouEXt3xydSDx5YNjghhEgvraFmTfi//zPac+aAdcxgQ3Q0eHvn6UQqLCqMXw79QtflXZMtbfBW7bcem0hBygsSy0LFIj3yTTIF0LdDX/p+3pc7R+/QfOJ7xsbNm43HgHNJgTkhRD725Zfw+uvGe6WMJ/LatTPmSs2blyfXx3vU2Ttn+d+m/1H2+7L0XdOXOw/vMKb5GOytky9tkJrt27cnmoSekCxULNIjXyVTAPO/nE+VNlXYs2QPn0z5BMqVg5dfhooVLR2aEEIktmsXeHgYw3kADg5QqFB8e/hwY26Ul1f8tlix6+PlAVHmKNacXEObRW2o+GNFJu+dTItyLdjSZwsnh55kdLPR9Hftj5WyiittkBqtNd988w0tWrSgZMmS2NvbJ9ovCxWLdNNaW+RVr149bSlBIUG6YLmCGlv0ht0b4ndERGg9frzWoaEWi00IkY/dvav1nDla37ljtBct0rpUKa0vXEj9PFfX5BZ0MbbnYlfvX9Vjt43Vjt87asagy35XVo/dNlZfvX81ybHX7l/TTec11YHBgale8+7du/q1117TgO7WrZsOCgrSixYt0s7OzloppZ2dnfWiRYuy6iOJXAzw1SnkNPlmAvqj9p/YzwsNXsCmoA3njpzDsaQjbNxodJmvWQOdOlksNiFEPnLlijHXqXRp2LcPXnjBeNq4Z0+IjAQrK6MuVB6VUmmDonZFeRD5gChzFK0rtOZdt3fpWKkj1lZPXtHnwIEDvPHGGwQEBPDNN9/w4YcfpjjMJ8SjZAJ6MupXrc+kWZOIuBlBgw4NMJvN0KYNHDkSn0iFhFg2SCFE3hQWZnwNCjKmGkydarTr1zfmb/boYbRtbPJ0IgVGaQNbU9J17sKiwhj2wjDOvH+Gjb030rlK5ydOpLTWzJo1i8aNGxMeHs4///zDsGHDJJESmSbfJlMAH/b8kM7vdyZwfyDth8Qsalm9uvH17FmjOvCKFZYLUAiR97zyCvTrZ7wvWhTmzoW+fY22UkZNqHzyS96szTR1akq0OTrRdhsrG04MPcE3rb/h+aeez9A9QkND6devHx4eHjRt2hR/f38aN26coWsK8ah8nUwBrPp+Fc5NnNk4cyNfz/86fkeJEsY/evXrWy44IUTuEhgIzZolforu55+hS5f4docO0LJlfPutt+D5jCUMuc2FuxcYvXU05X4oR69fe2GyMsWVMbA12fJO3XcoV7xchu9z+vRpXnjhBRYuXMjo0aP5448/KFmyZIavK8Sj8n0yZWVlxb51+7AvY4/nUE92HNph7ChWDHx8wMXFaE+ZkicfMxZCZCIvL2P9u9atjfIEYDxVFxERP7Q3fDi8847lYrSQh5EP8TnsQ8tfWlJ+Snm8tntRtURVlnVbxsmhJ+OG+kzKlObSBqlZsWIFbm5uBAYG8scffzBmzBhMeXzIVFhOvk+mAEoVL8WGdRvADG1ebcPtoNuJD7hwAUaMMNa4EkKIhG7cgMmTwdfXqO2ktTH3ckfMH2bvvw8bNsAjj9/nB1prfK/5MmTDEMp8V4beq3tz4e4FvF724uKwi/zZ+0+6V+9OueLl0lXa4FE+Pj64uLhgZWWFs7Mzbdu2pXv37lSvXh1/f3/atGmTRZ9QCEO+fZovOWNnjWWMxxjKNy/PmS1nsLJKkGueOmV0xZtMEBxsLN8ghMh/QkNh8WKoW9d4nTplLCDcsqWRQEVEgK2tUVAzdmJ5HpbS03hlC5flqQJPceTGEeyt7elWrRsDXAfQzKVZspXJA4MD6bmqJ8u6LUtXMvXo2nqx2rRpw7p167C1TTq5XYgnkdrTfJJMPaJwm8KEbAoBeyAMKAq0BNfWrvgP8of7941Hl3v2hNGjLRytECLLaW08iPLUU8Y8yocPoXhx+N//YNw4Y7+fH7z0UvxQHhjr450/b5Q8yMOGbBjCHP85SRYaBmhQtgEDXAfQs0ZPitoXzZL7y9p6IrtIaYR0aFCjASiMRAogCPgNSpwpYbQLFIC2baF5c8sEKITIelu2wMqVxnulYNQomDbNaBcoYPRGxVYXV8p4Ii8PVyBPidaajpU6Jnkaz1pZs7XvVv4d+C+D3AZlWSIFsraeyBmevPpZHnVqxSl4tLMuEk4uPwnfYtR9mTQpft/ixUY5hdq1szNMIURm8vc3XgMGGO0ffoBz56BbN6O9aROULRt/vLNz4vP37ImfcB4rIgJ27866mC3o9O3T+Bz2YfHRxZy9cxarmP+ZMWNrsmVgnYE0d2me5XGcPXsWGxsbIh793iNr64nsJT1Tj7h25Vqy269cvsIM3xncDk0wOT0sDDw98/xfn0LkOefPw08/GUN0YPxR9N57EB5utH/+2ZhQHsvZ2ahSnhJ//+QWc8lTC6hfD7nO5L2TqT+rPpV/qozXdi+cijoxp9Mcjg09hq115j6N9zg+Pj7UqVMHk8mEnZ1don2ytp7IbpJMPSKlv2asi1szeMNgSn9Xmg6LO7Do8CKCVaTxF2nsU34PHhhd+0IIy0iuzhMYT9zNnAn37hntv/82nrI7c8Zo/+9/cPUqxP5SdnQ0hvPygToz6qDGqiSvOjPqcD/8PvMPzqf1wtaU/b4swzcOx6zNfNf6O6783xW29NnCgDoDqFKiSoaexkuP4OBg+vbtS+/evXF1deXkyZPMmTMHZ2dnlFI4Ozszc+ZM3N3dsywGIR4lE9AfkdKTIXXq1GHG2hmsOLmCpUeXcvn+ZQpYF6BjpY70qtGLdhXaYN/lDerU2MlBh/tJrutaOmYCuxAi6wwZAjNmwNtvw8svQ4MGxkoGO3caE8TXrYNXXzWWcQkJSTx0l08lN4Hc2soa56LOXA2+SlhUGOWLl+fNGm/yZs03qVqyarLXedKn8dLDz8+PXr16ce7cOUaNGsXIkSOxTq3HUIhMJE/zpZOPjw+ffvYpVy9fxfE5R1q2aMmCBQt4++23mTVrFhrN7su7WXJkCSuOr+Bm6E2K2BWhq1V1bkWHsCn6VKJ/mGLnEEztkPcfkxbCIiIjjXlOnp7GXCV7e2MYfuJE+PhjY//581CpUr5ZqiWtAoMDKT+lPGFRYYm2P13gaXrV6IV7LXdeKPuCRdexM5vNTJ48mREjRvDMM8+waNEimjVrZrF4RP4kyVQmGDVqFOPGjWPkyJF4JZgjFWWOYsv5LSw5uoTVJ1dzPzxpr1QBkz3nh13I0q5vIfKdjz825jK9954xP8neHqKijKfqbG2NJVx8fPL8QsFPKiI6gr/O/8Wq46vwOeJDeLQxX8wKK9pVbMfqHquxMdlYOEq4ceMG/fr1448//qBz587Mnj2bp59+2tJhiXwotWRK+kfT6Msvv+T69euMGzeOMmXKMGTIEMDoDm/zfBvaPN+G6VHT+f3M73y8uB/nTcFGiQUNtR4Wxfx/w6Hws/Ddd8YF//0XSpUyVowXQiRlNhtznWLrNPXpY/Q6LV1qtA8eNHqcwJgjpVR8eYKICGNI7+bNPF/nKT0eRj5k47mNrDqxit9O/UZQeBBF7IrQoVIH1p9eT0R0BHbWdszuNDtHJFJ//fUXb731Fnfv3uXnn39m8ODBFu0hEyIlkkylkVKKadOmcePGDd577z2eeeYZXn/99UTH2Fvb07VYIxpNi6D8YAizASsN/9r+x3PPLqOtuTwDjjfm1cqvYtu7t1E9edky4+Q+faBhQ2POBxjLUTz3nLFGoBD5QUAAHD9u1HEDozDukSNw4oTRrlw58QMemzfHv/fyin8yL1Zsnad8XIXctbQrO/rv4Pczv7PqxCo2nN7Ag8gHPFXgKV6v+jqvV3udluVaYmdtx5ANQ5jhNyPLJ5CnRWRkJKNGjWLixIlUrVqVTZs2UbNmTYvGJERq5Gm+dLC2tmbJkiU0atSIN998k3/++SfpQV5elAnW9D8IVmYY7Aunp9kwIrQOh4qF0W1FN5797lmGedbj0JCuxjlaw7VrxqRYMP66btAgvuSC1uDuDuvXx7cvXjSGNFKS0lNNQmSl9Pzc+frCyJHxSdCPP0LnzvG9TQMGwGefxR/v6QlffJH8tfJZnadHNXJsFLdQcCxrK2vuhd2j5Dcl6bGyB9subqN3rd5sfmsz1z+6zpzX5tC+YnvsrI0nGEc1HUUTpybZUtbgUQnX1itbtixVqlRhwoQJeHh4sH//fkmkRM6ntbbIq169ejq3un37tq5WrZouUqSIPnToUOKdrq5ag75WCN20HzqwUEzFGVdXHRUdpf8484fuvqK7tvWy1YxB151RV//474/6dujt+GtERWm9Zo3WBw8a7Xv3tK5USevp0432jRvGNSdNMtp372o9ZIjW+/cb7chIrQcO1NrKytguRHZ5993EP3dms/HSWuvt27Vu1cr4+dVa6xkztLa11fryZaN97pzxMx8dnf1x53LX7l/T9l72mjEkepX+prR+//f39T8X/9FR0VGWDjNZixYt0g4ODhqjXHLc64MPPrB0aEIkAvjqFHIamYD+hC5fvkzjxo2Jjo5m9+7duLi4pOv826G3WXJ0CXP95+J/3R9bky2dq3TG95ov5++eT3J8otIK9+/D8uXQuDFUq2YMjbz4orFifefOsHFj/FBJgQLGX+1+fsYj4SVLZvCTC5GCy5ehYkWj8GWBAsZ6dv37w5o1xs/qjh3wwQewcCHUqGGscWdlFV/bSaSL1hr/6/5sOL2B9WfWs+/qvrh9JmXitcqvsaL7imQXFc5JZG09kVvI03xZ5NixYzRp0oRSpUqxa9cuSpQo8UTXOXj9IPP85+FzxIfbD28n2Z/m0gpms/HL6a23YMkSY86Ira1Rb2fjRiOhqlsXtm6Fb76B6dPByQmCg43zChZ8ovhFPvLggTG8XLQo3L0LI0ZA9+7QsiW8+abxcwfGz12vXsbP1bBhUKuWZePOIx5EPOCv83+x4cwGNpzZwLXgaygULzi+QDOnZkz+dzLh0eEUsC7A+Q/PW3zuU1pYWVmR3O8hpRTmR9c7FMKCMrzQsVKqrVLqlFLqrFJqRDL7nZVSW5RSh5VS25RSjhkNOjeoXr06v/32GwEBAXTo0IEHDx480XVcS7vyQ7sfuPp/V5nZcWaSvyTN2ky3at2S/QcnESsrY87KypXxE3UjImD7dqN3qkYNY1twsFHtuUgRoz13LhQqBLduGe1//zV6D2Lnroj8JXYuntbg7Q0bNhjthw+hcGFjbhMYvU9r1hj1mwIDYfXq+GtERBi9p199JYlUGqRWhfzivYv8tO8n2vm04+mJT9N5WWeWHVvGi8+9yILOC/jvctpXeAAAIABJREFUf/+x5+09jG81ngF1BmRLFfLM8PDhQ4YNG5biv2uytp7IVVIa/4t9ASbgHFAesAUOAdUeOWYF0DfmfQtg4eOum5vnTD1qzZo12srKSrdr105HRERk+Hrvrn83bk6VGqPi5j9UnFJRj9g8Qu+/ul+bY+ehJDn5XWMeSsIVwmxtU5875een9fjx8XNb3n9f60KF4tvffKN1jx7xx9+4oXVYWOof4to1rZs21TowMO0fXGSO9HzvV6zQ+rff4tuVK2vt4RHfLl1a6+HD49uTJmm9d2/S6zzJz52Ik/D/87Ev01iTLj6+eFy70o+V9P/9+X/67/N/64io5P+duXb/mm46r6kODM7Z/7/bt2+frlKligZ0mzZtdIECBRLNl3JwcNCLFi2ydJhCJEIqc6bSkkw1AjYmaH8GfPbIMceA52LeK+D+466bl5IprbWeOXOmBnTfvn1TTnTS6Nr9a9p+nDGZtMC4AvrQ9UN62v5p+pVfXtGmsSbNGLTzJGc9/M/heselHTranGDCbswE+CQvV9e0BxARofX58/Htr77SukuX+HaXLlpXrRrfXrFC6w0bEl/j0YnIIvsk/N7v36/1xo3x+3r10rp79/i2m5vWrVvHt8eP13r58vh2eHja7pkZP3f/3959x0dVpQ0c/52ZtAkt1BBKhrorCAKiiCKsoC7FRUR0FXglKoKCuqvivqBhlWLsIvpaUbFscMWlLSIL0talKE1ADbgSSmihtwBpk3neP2ZSJsmEhEkyM8nz/Xzmk9vm3ufMnZk8c86551ZjGw5skNApoUU6kF//0fXy+nevy6/Hf/V3iOUiKytLnn32WbFardKsWTNZvny5iLg6odvtdjHGiN1u10RKBSRfk6k7gA8LzN8DvFVom8+BP7unb3f/uqhfzL5GA5uATbGxsZX3ClSSyZMnCyDjx4/3eV9jFo0Ry2SLjF3kmYycuHBCPt7ysfzh8z/k/ZJt/GpjGbNojCzftVw6vdupyBcyk5DO75XjP7XFi0Vmzcqf79xZZMCA/PnBg0VCQlxvL5tNZPp0kW+/zV+/c6fIyZPlF09h/q4Vq8jj790rsmZN/vyCBSIJCfnzjz3mSqRyX/v+/T0T3+ef99w+NfXitYyq3J24cELmJM2Rh756SNq82abI5zVkSoiM/OdIf4dZrrZv3y5du3YVQEaMGCGnTp3yd0hKlUllJFNNgHnAFuAN4AAQVdJ+q1rNlIiI0+mUhx56SACpW7euT7+ySlNdfybjjHz+4+dyx5d3SGRCpDAJCZ8aLpZJFo8v5rCpYUWSsnKVnu6ZOLRtK2K15jf1RESIPPhg/vq6dUUefjh//re/ddWI5Bo1SmThQte00yny+eciv/ySP3/0qGv4B2/8XStW0vGzs0UOH86//H/7dpGZM/Pn580TufPO/CbWKVNcr1euxx8XqVEjf/6RR0SaNMmf79EjP5kKCxMZPlxkz55yLZ4qqvN7nUv8EXMh64J8k/yNjF82Xrq+3zWv+b7W87Vk4OcD5Y3v35BVu1d51EgHelNdaeXk5Mj06dMlIiJCGjRoIHPnzvV3SEpdkgpv5iu0fU3gwMX2WxWTKRGRzz77TKxWa6W3/5/POi/zts+T27+4vcgXeuiUUFm9d7XPzY+lcuiQK3kq2NQTEeFKGnLNni2yfr1r2ukUGTnStUzE1azUvLnIq6+65tPSXPt46SXX/OnTrvnXXnPNnzjhqhnL/YLesSM/kbPZRH76SWTqVJGkpPz4Jk/OT8727xd59lmRX93NKHv2iPz1r64xj0REkpNF4uNdNUIiruc99ZTIvn2u+Z9/Fhk/3rVfEZFlyzxr5T77zNWUduSIa/2bb7rW5Y61NG2aaz63pu7dd101SefPu+a//lrkL39xjT2We/x//zs/2So4JlNxr73Npv3WKkFxfZ5Cp4TKNR9cI30+7SPhU8PzlvX6uJdM+fcUWbtvbZG+T95qpINVSkqK9O7dWwAZOHCgpOp7UQUxX5OpEGA30LJAB/TLC23TALC4pxOAKRfbb1VNpux2e5HB5wCx2+2VFsPohaMlZEpIkaTK/rpdRi8cLXOS5sip9AqqYi/vjsgOhytByk1Gzp1zJSRbtrjmDx8WGThQZMkS1/zQoZ7Hvftu13RusvbDD675BQtc899/75rP7fP17bcixoisWOGaX77clZytXu2aX7xYJDRUZMMG1/yCBSLh4fkDrP7+957HHzhQpF8/kQMHXOt//FHkrbdEzp51zZ844UrgHOUwoKJ2Aveb/af35yVMhR+d3u0k45aOk8W/Lpa0zLQS9xMsHcgvxul0yieffCK1a9eWmjVrykcffVQ5P+aUqkAlJVOlGmfKGDMAmI7ryr6ZIpJgjJni3vFCY8wdwAvuxOE/wMMiklnSPqvCOFPFCYQxU1LTUmn1ZisyHBnYQmysilvFlsNbWLprKSt2ryAtKw2rsXJNs2vo27ovfVv35aomV2G1WEu8x1feoKEl6dLFdQPaIjvoDFtK8XxfpKZCq1aQkZG/zGaDX3+FJk1cQ0fkphnGuB650xV5/N27K+dmu/587auZ81nnWX9wPWv2rWHNvjV8f+B70rLS8tZbjIUbW95I4u2JNKrRyI+RVo5Zs2YRHx/Pvn37aNq0KdHR0WzevJlevXrxySef0FJv6K6qAB20sxJ5G823Xr16nDhRdEDOipJ709KHuj7kMdhndk423x/4nqW7lrJ011I2H9qMINSz1eOmVjdx9PxR1u1fR1ZO/n3OSj1oqL+NHQsffeR5j7awMHjggcq52a2/j68u2cV+RKSmpbJ2/1rW7lvLmv1r2JK6hRzJwWDoGN2RHs17cHnDyxn3zbigGzTTV7NmzWL06NFcuHDBY/nQoUNJTEzEYgnsEdiVKq2SkqmQyg6mqktISCjyxWKxWDh58iTTp0/nscceq5Q4/trrryQdSypy09JQayg97T3pae/Jc32e4/iF4yzbtYylu5byza5vSD2XWmRfFmPxy81Py8zfN7v19/HVJbu22bVsP7bd40eE1VhJz06n9Zut827xZAux0a1pNyZcP4EezXtwbfNriYqIyntO0rEk3t/8flAMmlle4uPjiyRSAOvWrdNESlUbWjNVAQpWecfGxjJ58mS++uor5s6dS0JCAk8//bS/QyyWiPDz0Z95cNGDfH/ge4T890az2s3o0byH6xHbgyuiryDEorm4Cm4iwsG0gyzbtYzRi0bjcDo81te31aeXvRfXx15Pj+Y96BLThTBrmNf9paalcvfcu5l9x+xqkUyJiNeESW8Ho6oabeYLAA6Hg/vuu4/ExETi4+OZOnUqprz66pSzgn2uwq3hTOw1kZ+P/sza/Ws5cPYAADXDanJN02vykqvuzbpTO7y2732ulLoEpXnf5SZOmw9tZnPqZjYd2sTm1M0cPX+0yPOsxspdl99F4u2JAfs59bf9+/czZswYvs693VAheqNiVdVoM18ACAkJ4dNPP8Vms5GQkMD58+eZNm1aQH5Rx9SK4b7O9/H+5vcZ2WUkE3tNzFu378w+1u5b6+o/sn8tz61+Dqc4sRgLHRt1xOF0EGIJ8fiFH2YN47pm1/mjKKqaKK6ZLswaRuMajXlm1TN5yVNu4mQxFi5veDkD2g6ga0xXusZ0pVHNRnR4pwMZjgzCrGG81ve1gPx8+pvT6eT9999n/Pjx5OTkMHz4cObPn+/R1BcZGUlCQoIfo1SqcmkyVYksFgvvv/8+NpuN6dOnk56ezjvvvBOQ/Qq89bmKrRNLbMdYhnYcCsDZzLOsP7A+L7lat39dkaYSp9PJbxv8lh+P/MhlDS4rsZlEqbK6kH2BQZcN4sMfPvRYnpWTxZJdS/hm9zdFEqdOjTsRGRpZZF+5PyKqU5+nsvjvf//LqFGjWL16NTfddBMzZsygZcuWRbo2JCQkMHz4cH+Hq1Sl0WY+PxAR4uPjeeGFF7jnnnuYOXMmISFVI691OB0MmzOMeb/My7vayWIs5EgOAKGWUNo3bE/nxp3pFN3J9bdxJ+rZ6gGla65RVU9pzvuF7Av8cvwXko4msf3YdpKOuf7uPrXbo38fgAUL18Vex8s3vew1cSpOdevzVFrZ2dm88sorTJkyhcjISKZNm0ZcXJzW3KlqRZv5Aowxhueff54aNWowceJE0tPTmTVrFmFhwV9jE2IJ4Y3+b/DVzq/IceQQERLBr4/+ytnMs2w7vI2th7ey7cg2lu5ayqfbPs17XvPazenUuBMWY9FmwmqouGa6EBOCwXDr328l6VgSe07tyUuaQi2h/Kb+b+japCsjOo2gfcP2NIpsRN9ZfV19/ULC+ced/yhzQhRTK4Zv7/22XMsW7DZt2sQDDzzAtm3buPPOO3nzzTdpXBnjpikVRDSZ8qP4+HhsNhvjxo0jIyODf/zjH0RERPg7LJ8V7HN1X+f7aFa7GQDtG7bPax4EOHLuCNuO5CdYWw9v5Zdjv+DE8wogh9NBhiODN9e/Sdt6bWlTrw0toloQag312E5rtfynLK+9U5wcPneYPaf2sPvUbnaf2s3xC8dx5Hg2DzvEwU9HfiLbmc1VTa4irlMclze8nPYN29OmXpsi5x+0ma48XbhwgUmTJvHaa68RHR3N/Pnzue222/wdllIBSZMpP3viiSew2WyMHTuWgQMHsmDBAmrUqOHvsHzmrc9VQdE1o/l9zd/z+9a/z1uWnp3OPfPvYcEvC8iRHCxYqBtRly+3f8nMrTPztrMaK/Yoe15y1aZeG2JqxpBkSSLbmZ23XVlqtap7MuZL+YurWQq1hNK4RmNe/+51V9J0ejd7Tu1hz+k9ZDjyR4k3GJrWbkp0zWiOnDuCEychlhDuaHcHnw3+rNikyZvSvO9UUYX7PA0fPpzZs2eza9cuRo0axcsvv0xUVNTFd6RUNaXJVAAYM2YMkZGR3H///fTr14+vv/6a2rVr+zssn1xqc4kt1Mb/9f8/vt75NTmOHMJDwvl57M9E14jm6PmjJJ9MJvlkMjtP7syb/u7Ad5zNPFvs/nKcOaQ70nl+9fM0qdWEmJoxxNSKoUmtJtS31ffo8+HtirBgScZ8PX5pyp+Wmcbhc4dJPZdKalpq3vTxC8eLXHiQ7cxmya4lLNm1hNrhtWlVtxXtGrZjQNsBtKrbilZ1W9EyqiX2KDsRIREeQ3KEWkJ5vd/rZUqkQJvpLkXhEcxTUlJ4/vnnadSoEStXrqR3795+jlCpwKfJVICIi4vDZrMxfPhwOnfuTHZ2NgcPHqyWV8YUbibMba6JrhlNdM1oesT28NheRDiRfoKdJ3YyceVEvk35Nq/ze62wWsz/ZT6nM04XOU6oJZTGNRu7kqxaMdQOq41TPJsYDYZ7O9/LqfRT1A6vjdVi9Rq3v5Oxsh4/05HJmcwznMk4w5nMM9xgv6HIFXE5zhw2pW6izZttOHzuMOezzxfZT+7rWN9Wn+MXjiMIVmPlxpY3knBjAi2jWlLPVu+inZW9nXdVsbyNYB4eHq6JlFKlpMlUAPnjH//I+vXrmTZtWt6ylJQURo8eDVCtEqqyNNcYY2gQ2YAGkQ1IvD2RVm+2yuv8vuORHTSu2Zj07PS82pRDaYfyp88dIjUtlZ0ndpJ6LrVI7UpmTibdPuyWN18zrCZ1wutQO7w2dSLqUCe8Tt5fi7EUScYQVzI0f8d8Qq2hhFhCCLW4/xaYD7WGclmDy0g66tlMGWoJ5Tf1fsOGgxvIzskmKyeLbKf7b6H5FnVaFLnJdo4zh5QzKdz02U0eidOZjDNk5pR4L/K88kaGRtKqaStiasbQuGbj/L+1XH/r2ephMRaPmqUwaxifDv60zAmRNtNVvn379hW7/MCBA5UciVLBS4dGCDDebpSsowmXnrebPJfGnlN7aPd2OzJzMgmzhvH2gLexGEteEnI282x+QuJOSs5mns2bTnekV1CpLk2N0BrYo+weSZ/HdEQdoiKi8qazHFn84e9/uOSb9fry2qvKlZaWxpQpU3j11VeLXa/fOUp50qERgoi3X4nelquifKndaFm3Jfd3uZ/3N7/PA10e4IErHyjT81NOp3DZ25eR4cggwhrBsnuWUddWl2xnNg6ng+ycbI9ph9PhMf/hDx/mNVOGWELo06IPY64eQ6gllDBrGKFW918v8yfST3D1B1eT4cjAFmIj+U/JZa4dyi3/pTS1ac1S4BMRvvzyS5544gkOHTrE7373OzZs2EB6ev4PAR3BXKkyEhG/PLp27SqqKLvdLkCRR926dcXpdPo7vGrh0NlD0uvjXpKalnpJzx+zaIxYJltk7KKxl3TsiOcihEmI7TnbJcXgy/FzY/Cl/CpwJSUlSZ8+fQSQLl26yHfffSciIomJiWK328UYI3a7XRITE/0cqVKBB9gkXnIaTaYCTGJiokRGRnokUlarVQCJi4uT9PR0f4eoLsKfyVh5HF9VPWfPnpUnn3xSQkJCJCoqSt555x1xOBz+DkupoFJSMqV9pgJQ4TFfnnvuOXbt2sWkSZPo3r078+bNIyYmxt9hqgqitzRR5UVEmD17NuPGjePQoUOMHDmSF154gYYNG/o7NKWCTkl9pjSZCiJz585lxIgR1K1bl/nz53P11Vf7OySlVIAo/CNs7NixLFmyhFWrVnHllVfy9ttv0717d3+HqVTQ0mSqCtm2bRuDBg3iyJEjfPTRRwwbNszfISml/KzwwJu5bDYbr732GqNHj8Zq9T5GmlLq4kpKpiyVHYzyTadOndi4cSPdunVj+PDhTJgwgZycHH+HpZTyI28Db9avX58xY8ZoIqVUBdNkKgg1bNiQZcuW8eCDD/LSSy8xaNAgzp4t/nYqSqmqz9vQKQcPHqzkSJSqnjSZClJhYWG89957vPPOOyxZsoTu3buTnJzs77CUUpVo586dDBkypMjI97liY2MrOSKlqidNpoLcmDFjWLZsGUePHqVbt24sX77c3yEppSrYsWPHePTRR2nfvj1Lly5lyJAh2Gw2j2104E2lKo8mU1VA79692bBhA02bNqVfv36MGDECu92OxWKhRYsWzJo1y98hKqXKQXp6Oi+++CJt2rTh3XffZeTIkSQnJzNnzhw++OAD7HY7xhjsdjszZsyoVvfzVMqf9Gq+KiQtLY0bbriBH374wWN5ZGSkfrEqFcScTieJiYlMnDiR/fv3M3DgQF566SXatWvn79CUqjZ8vprPGNPPGPNfY0yyMWZCMetjjTGrjDFbjDE/GmMG+Bq0KrtatWpx/PjxIssvXLhAfHy8HyJSSvlq+fLldO3albi4OBo1asSqVatYuHChJlJKBZCLJlPGGCvwNtAfaA8MNca0L7TZROBLEekC3A28U96BqtLZv39/scv1RslKBbZZs2bRokWLvOb5l156iQEDBnDzzTdz6tQpZs2axYYNG7jhhhv8HapSqpCQUmzTDUgWkd0AxpgvgEHA9gLbCFDbPV0HOFSeQarSi42NJSUlpcjy8PBw9u/fT/Pmzf0QlVKqJIUH3UxJSWHChAnYbDZeeeUVHnnkESIiIvwcpVLKm9I08zUFClZ3HHAvK2gS8D/GmAPAYuDR4nZkjBltjNlkjNl07NixSwhXXUxCQgKRkZEey8LCwnA6nXTs2JG//e1vXi+jVkr5R0mDbj755JOaSCkV4Mrrar6hwCci0gwYAPzNGFNk3yIyQ0SuEpGr9EabFWP48OHMmDHD46qemTNnsn37djp06MCIESO444470GRWqcBw+vTpYmuTQQfdVCpYlCaZOggUbBtq5l5W0EjgSwAR+Q6IABqUR4Cq7IYPH87evXtxOp3s3buX4cOH07p1a7799lteeuklFi1aRMeOHVm0aJG/Q1Wq2jp9+jSTJ0+mRYsWXrfRQTeVCg6lSaY2Am2NMS2NMWG4OpgvLLTNPuBGAGNMO1zJlFZ9BBir1cr//u//snHjRqKjoxk4cCCjRo0iLS3N36EpVW0UTKImTZpEnz59im2e10E3lQoeF02mRMQBPAIsBXbgumovyRgzxRhzq3uzccAoY8w24O/AvaIdcwLWFVdcwYYNG3jqqaeYOXMmV1xxBf/5z3/8HZZSVVpxSdSWLVuYN28eTz/9dJHmeR0bTqngoYN2VnNr164lLi6O3bt3M27cOKZOnaqdXZUqR6dPn+aNN97g9ddf58yZMwwePJhnnnmGzp07+zs0pVQZ+Dxop6q6evTowdatW3nwwQd59dVXueqqq0hISPAY70ZvR6PUxRUeJ+qDDz7wWhOliZRSVYvWTKk8//rXvxg6dChnzpzxWK63o1GqZIXHiSpIa6KUqhpKqpnSZEp5aN68OQcOHCiy3G63s3fv3soPSKkg0KJFi2KHN4iJieHQIR3DWKmqQJv5VKl5G9dGb0ejVPGSkpK8jhN1+PDhSo5GKeUPmkwpD97GtRER/vznP3PixIlKjkipwCMifPPNN/Tr148OHTpgjCl2Ox0nSqnqQZMp5aG48W5sNhs33ngjb731Fm3btmX69OlkZWX5KUKl/CcjI4OZM2fSsWNH+vbty7Zt20hISODdd9/VcaKUqsY0mVIeirsdzQcffMDy5cvZtm0bV199NY8//jgdOnTgn//8p97nT1ULx44dY8qUKdjtdkaOHElISAiffvope/fu5emnn+bBBx/UcaKUqsa0A7oqExFhyZIljBs3jh07dtC7d2+mTZumVyqpoDdr1izi4+PZt28fsbGxJCQk0KVLF6ZPn85nn31GZmYmt9xyC0888QS9e/f22rSnlKqa9Go+Ve4cDgczZszgmWee4eTJk9x3330899xzxMTE+Ds0pcqsuKENLBYLTqeTiIgI4uLieOyxx7jsssv8GKVSyp80mVIV5vTp0yQkJPDGG28QFhbGhAkTaNq0KZMnT/b4ha/NHSqQeRvaoE6dOiQnJ9Oggd63XanqTpMpVeF27drF+PHjmTt3LsYYj75UOuinClQiwtq1a+nZs2ex640xOJ3OSo5KKRWIdJwpVeFat27NnDlziI6OLtIp/cKFC8THx/spMqWKOnLkCC+//DLt2rWjZ8+eOrSBUsonmkypcnX06NFil6ekpHhdp1RlcDgcLFq0iMGDB9OsWTPGjx9PgwYNmDlzJh9++KEObaCUumSaTKlyVdIv+djYWEaNGsWOHTsqMSJVXRS+0XDuDbqTk5N5+umniY2NZeDAgaxbt47HH3+cHTt2sGbNGu677z7uv/9+HdpAKXXJtM+UKlfFXRUVGRnJlClTSE5O5pNPPiEjI4MBAwYwbtw4vcRclYvi3ndhYWG0atWKX375BYvFQv/+/XnggQe45ZZbCA0N9WO0SqlgpB3QVaUqbrye3F/4x48f59133+Wtt97i6NGjdO7cmSeeeIK77rqLsLAwP0eugpW3q/FCQkKYPHkycXFxNG3a1A+RKaWqCk2mVMDJyMhg1qxZTJs2je3bt9OkSRP+9Kc/MXr0aBYvXuw1GVOqIIfDwerVq+nTp0+x6/VqPKVUedFkSgWs3BHVX3vtNVasWEFYWBhOpxOHw5G3jQ6toArKyspi5cqVzJ07lwULFnD8+PEiw3Hkstvt7N27t/KDVEpVOTo0ggpYxhj69+/P8uXL2bJlC6GhoR6JFOjQCtWFtw7k4KrJXLhwIXFxcURHR9O/f3+++OILbrrpJubMmaNX4yml/CrE3wEolatz584eHYgLSklJ4aOPPmLw4MHUq1evkiNTFa1wB/KUlBRGjRrF+vXrOXr0KF9//TXnzp0jKiqKQYMGMWTIEG6++WYiIiLy9hEeHq7Nw0opv9BmPhVQSupI7HA4CAkJ4eabb+buu+9m0KBB1KlTxw9RqvLm7bwDNGzYkNtuu40hQ4bQu3dvvVBBKeUX2syngkZCQkKxzTWffPIJGzdu5PHHHycpKYm4uDgaNWrEbbfdxueff865c+eAkpuKVOA5e/YsixYt8ppIGWM4dOgQM2bMoG/fvppIKaUCktZMqYBT0tAK4Oq0vn79embPns2XX37JoUOHiIiIoGPHjmzbto2srKy8bbXzemBJT0/nu+++Y8WKFaxcuZKNGzeSk5PjdXvtQK6UChR6NZ+qspxOJ2vWrGH27Nm89957xV4GHxsb67XmQ5Wf4pLgP/7xj2zatImVK1eyYsUK1q1bR2ZmJlarlW7dutGnTx/69OnDvn37ePjhh4sM9qqJsFIqUGgypaoFi8VS7OXxALfeeis9e/akZ8+eXHnllToCdjkrbgRyq9VKSEgImZmZgOsCg9zkqVevXtSqVavIPrQDuVIqUPmcTBlj+gFvAFbgQxF5sdD614He7tlIoJGIRJW0T02mVHnz1om5Ro0axMTEkJycDIDNZqN79+55yVX37t2pWbMmoP/Qy+Ls2bP8+OOPbN26lQkTJnD+/Pki29SsWZOPP/6YG264gQYNGvghSqWUKh8+JVPGGCvwK3AzcADYCAwVke1etn8U6CIi95e0X02mVHnzdl/A3Kaiw4cPs2bNGlavXs3q1avZtm0bTqcTq9XKlVdeSf369Vm1alVeTUrh51cHxSWTw4YNIzU1lS1btrB169a8R25yWhIdgVwpVVWUlEwhIiU+gGuBpQXmnwKeKmH7dcDNF9tv165dRanylpiYKHa7XYwxYrfbJTEx0eu2Z86ckSVLlkh8fLz06tVLgGIfdevWlcWLF8uvv/4qWVlZ5Xb8iuDL8RMTE8Vms3mU3WKxSK1atTyWtW7dWoYMGSJTp06VRYsWyYEDByQ2NrbY185ut1dcYZVSqhIBm8RLTlOamqk7gH4i8oB7/h7gGhF5pJht7cD3QDMR8X6JDlozpQJPSX2uclmtVux2O23atKFt27a0adMm77F+/XrGjh3rUydqX5oZL1Yzl5GRwf79+70+kpKSii1/jRo1eOGFF+jSpQtXXHEFtWvXLvOxlVIq2PnazFeWZGo8rkTqUS/7Gg2MBoiNje2qV1ipQOKtz1XsBomyAAAJX0lEQVSzZs344osvSE5O9njs3LmTM2fOXHS/UVFRTJ06lVq1apX4mD17tteEZNiwYWRlZZGenk5GRgbp6elFpocNG8axY8eKHD80NJSoqKhi1zVo0IDmzZvTvHlzFi5cWGz8pW2q0/5mSqmqzNdk6lpgkoj0dc8/BSAiLxSz7RbgYRFZd7GgtGZKBZqy1q6ICCdPnmTnzp0kJydzzz33VEhcxpi8412q0aNH5yVNuY9mzZphs9nytvGWTOpYT0op5XsyFYKrA/qNwEFcHdCHiUhSoe0uA5YALaUU3/qaTKlA5EvtirdkpHnz5mzatIm0tLQSH88++6zXfU+cOJGIiAhsNhs2m63Y6bvuuovDhw8XeW5pkyFtqlNKKe986oDuzosG4EqodgHx7mVTgFsLbDMJeLE0+xPtgK6qoMTERImMjPTogB0ZGVnqTuB2u92nTty+Hj93H/7sQK+UUoEKXzqgVxStmVJVUUV2IK/o4yullPJOR0BXKkhoMqSUUoFJkymllFJKKR+UlExZKjsYpZRSSqmqRJMppZRSSikfaDKllFJKKeUDTaaUUkoppXygyZRSSimllA/8djWfMeYYUBk352sAHK+E4wQiLXv1VZ3LX53LDtW7/Fr26qsyym8XkYbFrfBbMlVZjDGbvF3KWNVp2atn2aF6l786lx2qd/m17NWz7OD/8mszn1JKKaWUDzSZUkoppZTyQXVIpmb4OwA/0rJXX9W5/NW57FC9y69lr778Wv4q32dKKaWUUqoiVYeaKaWUUkqpCqPJlFJKKaWUD4I2mTLG9DPG/NcYk2yMmVDM+nBjzGz3+vXGmBYF1j3lXv5fY0zfyoy7PJSi7E8YY7YbY340xqwwxtgLrMsxxmx1PxZWbuTloxTlv9cYc6xAOR8osC7OGLPT/Yir3Mh9V4qyv16g3L8aY04XWBfU594YM9MYc9QY87OX9cYY86b7tfnRGHNlgXVBfd6hVOUf7i73T8aYdcaYTgXW7XUv32qM2VR5UZePUpT9BmPMmQLv72cKrCvxMxPoSlH2vxQo98/uz3k997qgPu8AxpjmxphV7v9pScaYPxezjf8/+yISdA/ACuwCWgFhwDagfaFtxgLvuafvBma7p9u7tw8HWrr3Y/V3mcq57L2BSPf0mNyyu+fP+bsMlVD+e4G3inluPWC3+29d93Rdf5epPMteaPtHgZlV6Nz3Aq4EfvayfgDwL8AA3YH1VeG8l6H81+WWC+ifW373/F6ggb/LUIFlvwFYVMzyMn1mAvFxsbIX2nYgsLKqnHd3GWKAK93TtYBfi/nO9/tnP1hrproBySKyW0SygC+AQYW2GQR86p6eA9xojDHu5V+ISKaI7AGS3fsLFhctu4isEpEL7tnvgWaVHGNFKs2596YvsExETorIKWAZ0K+C4qwIZS37UODvlRJZJRCR/wAnS9hkEPCZuHwPRBljYgj+8w5cvPwiss5dPqhin/tSnHtvfPm+CAhlLHuV+swDiEiqiPzgnk4DdgBNC23m989+sCZTTYH9BeYPUPTFzdtGRBzAGaB+KZ8byMoa/0hcGXuuCGPMJmPM98aY2yoiwApW2vIPcVf3zjHGNC/jcwNVqeN3N+22BFYWWBzs5/5ivL0+wX7eL0Xhz70A3xhjNhtjRvsppop2rTFmmzHmX8aYy93Lqs25N8ZE4koU5hZYXKXOu3F11+kCrC+0yu+f/ZCK2KkKDMaY/wGuAn5XYLFdRA4aY1oBK40xP4nILv9EWGG+Av4uIpnGmAdx1VD28XNMle1uYI6I5BRYVh3OfbVnjOmNK5m6vsDi693nvhGwzBjzi7vGo6r4Adf7+5wxZgCwAGjr55gq20BgrYgUrMWqMufdGFMTV6L4mIic9Xc8hQVrzdRBoHmB+WbuZcVuY4wJAeoAJ0r53EBWqviNMTcB8cCtIpKZu1xEDrr/7gb+jSvLDyYXLb+InChQ5g+BrqV9boArS/x3U6i6vwqc+4vx9voE+3kvNWPMFbje84NE5ETu8gLn/igwn+Dq2nBRInJWRM65pxcDocaYBlSjc0/Jn/mgPu/GmFBcidQsEZlXzCb+/+z7o0OZrw9cNWq7cTVj5HYqvLzQNg/j2QH9S/f05Xh2QN9NcHVAL03Zu+DqdNm20PK6QLh7ugGwk+DrjFma8scUmB4MfO+ergfscb8Odd3T9fxdpvIsu3u7y3B1PDVV6dy7Y2+B907It+DZCXVDVTjvZSh/LK4+oNcVWl4DqFVgeh3Qz99lKeeyN859v+NKGPa53wel+swE+qOksrvX18HVr6pGFTzvBvgMmF7CNn7/7AdlM5+IOIwxjwBLcV2tMVNEkowxU4BNIrIQ+Aj4mzEmGdeb7G73c5OMMV8C2wEH8LB4NoUEtFKW/RWgJvAPV5979onIrUA74H1jjBNXreSLIrLdLwW5RKUs/5+MMbfiOr8ncV3dh4icNMZMBTa6dzdFPKvEA1opyw6u9/oX4v42cQv6c2+M+Tuuq7YaGGMOAM8CoQAi8h6wGNdVPcnABeA+97qgPu+5SlH+Z3D1C33H/bl3iMhVQDQw370sBPhcRJZUegF8UIqy3wGMMcY4gHTgbvf7v9jPjB+KcMlKUXZw/Wj8RkTOF3hq0J93tx7APcBPxpit7mVP4/rxEDCffb2djFJKKaWUD4K1z5RSSimlVEDQZEoppZRSygeaTCmllFJK+UCTKaWUUkopH2gypZRSSinlA02mlFIBzxgTZYwZ655uYoyZ4++YlFIqlw6NoJQKeO57ci0SkQ5+DkUppYoIykE7lVLVzotAa/egfTuBdiLSwRhzL3AbrhGe2wKv4hrp+h4gExjgHrivNfA20BDXoH6jROSXyi+GUqoq0mY+pVQwmADsEpHOwF8KresA3A5cDSQAF0SkC/AdMMK9zQzgURHpCjwJvFMpUSulqgWtmVJKBbtVIpIGpBljzgBfuZf/BFzhvtv8deTfXglc9+ZUSqlyocmUUirYZRaYdhaYd+L6jrMAp921WkopVe60mU8pFQzSgFqX8kQROQvsMcbcCWBcOpVncEqp6k2TKaVUwBORE8BaY8zPwCuXsIvhwEhjzDYgCRhUnvEppao3HRpBKaWUUsoHWjOllFJKKeUDTaaUUkoppXygyZRSSimllA80mVJKKaWU8oEmU0oppZRSPtBkSimllFLKB5pMKaWUUkr54P8BeBx8CtebRaAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "ClSpSS27P7C7" }, "source": [ "## Convergent \n", "The Abysmal-Butler method does satisfy the Lipschitz condition:\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)=\\frac{4}{2}f(t,w_i)-\\frac{3}{2}f(t-h,w_{i-1}))-(\\frac{4}{2}f(t,\\hat{w}_{i})-\\frac{3}{2}f(t-h,\\hat{w}_{i-1})))\\end{equation}\n", "\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)=\\frac{4}{2}(f(t,w_i)-f(t,\\hat{w}_i))-\\frac{3}{2}(f(t-h,w_{i-1}))-f(t-h,\\hat{w}_{i-1})))\\end{equation}\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)\\leq\\frac{4}{2}L|w_i-\\hat{w_i}|+\\frac{3}{2}L|w-\\hat{w}|\\leq \\frac{7}{2} L|w_i-\\hat{w_i}|\\end{equation}\n", "This means it is internally convergent,\n", "\\begin{equation}|w_i-\\hat{w_i}|\\rightarrow 0\\end{equation}\n", "as $h \\rightarrow 0$,\n", "but as it is not consistent it will never converge to the exact solution\n", "\\begin{equation} |y_i-w_i| \\not= 0.\\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "nAWe9PM5P7C7", "outputId": "5e51f857-a073-46f3-827e-5461a8149d5f", "colab": { "base_uri": "https://localhost:8080/", "height": 203 } }, "source": [ "d = {'time': t[0:5], 'Abysmal Butler': w[0:5],'Abysmal Butler w2 N*2': w2[0:10:2],\n", " 'w-w2':np.abs(w[0:5]-w2[0:10:2])}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": 9, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
timeAbysmal ButlerAbysmal Butler w2 N*2w-w2
00.0001.0000001.0000000.000000
10.1250.9252230.9252230.000000
20.2500.9126670.8925950.020072
30.3750.8970420.8691120.027930
40.5000.8907820.8539750.036807
\n", "
" ], "text/plain": [ " time Abysmal Butler Abysmal Butler w2 N*2 w-w2\n", "0 0.000 1.000000 1.000000 0.000000\n", "1 0.125 0.925223 0.925223 0.000000\n", "2 0.250 0.912667 0.892595 0.020072\n", "3 0.375 0.897042 0.869112 0.027930\n", "4 0.500 0.890782 0.853975 0.036807" ] }, "metadata": {}, "execution_count": 9 } ] }, { "cell_type": "code", "metadata": { "id": "vae0CQb0P7C9" }, "source": [ "" ], "execution_count": 9, "outputs": [] } ] }